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ABSTBACT 

This thesis treats the problem of determining the eleva- 
tion angle needed to initialize an underwater sound ray 
tracing algorithm used to locate the position af a target 
vehicle. At regularly spaced time intervals the vehicle 
pings a synchronized sound signal which is received by a 
(short base line) sonar array containing four hydrophones 
positioned at four of the corners of a cube. Tie wavefront 
direction angles are determined from the arrival times at 
the four hydrophones. 

Current methods for using such time data ti produce an 
apparent position suitable for ray tracing are reviewed. 
Then four new methods are developed and documented mathemat- 
ically. All methods are compared under a simulated environ- 
ment of a sound speed profile which is linear with depth. 
One of the new methods is judged to be an impri vement over 
current methods in this idealized environment. Finally the 
improved method is used to estimate the variability in the 
time data from a real hydrophonic tracking problem. 
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I- INTEQDOCTION and background 

A. BACKGROUND 

Suppose that an underwater acoustical source emits a 
signal at a known time. Suppose that the times of arrival 
of that signal at the hydrophones on a three dimensional 
sensing array are also known. Finally suppose that the 
speed of sound in water is modelled as being homogeneous 
over time and horizontal displacement, varying only with 
depth. Then the angle of elevation (A), and the time (T) of 
the arrival of the signal at the acoustic center of the 
hydrophone array can be determined. If the relationship 
between the speed of sound and the depth under water is 
known exactly, then the angle A and the time T can be used 
to trace the signal trajectory back over its ray path 
(called ray tracing) to determine the original position of 
the acoustical source [Ref. 1]. However, full realization 
of this method in actual hydrophonic tracking Ls prevented 
by two primary sources of inaccuracies in the process. The 
first and probably greatest problem is that ta e speed of 
sound profile can be approximated at only a ferf locations, 
usually at great cost to achieve even moderate accuracy. In 
addition the profile is certain to flucuate over time and 
location. The second problem, confounding the first, is 
that there may be innacuracies , of unknown size, in the time 
data values recorded by the sensing array. 

B. PURPOSE 

As noted in [Ref. 1] the procedure of determining a 
sound source position by ray tracing is very sensitive to 
even small errors in the angle of elevation or speed of 
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sound at the sensing array. Of course ray tracing is also 
highly dependent on the accurate determination of the the 
transit time from source to array. The purpose of this 
study is to develop an appropriate model of the timing data 
observed in the hydrophonic tracking problem. Tie objective 
of the desired model is to estimate the error in the time 
values, and produce improved estimates of the iaitial angle 
and ray transit time, so as to reduce the ef f s ct of these 
errors on the ray tracing procedure. 

In pursuit of these objectives this thesis first reviews 
the currently used models, which are called the NAVY and 
NAVY_A methods. Then four alternative models are developed, 
called the L.S., L.S.C., M. L. P. and M.L.S. methois. Finally 
the performance of all six models are evaluated and compared 
using simulation studies. The comparisons are made under 
the idealized conditions of a known linear spe 2 d of sound 
versus depth relationship. 

C. TEACKIHG RANGE CCNFIGOR ATION 

The tracking range which supplied data for this study 
consists of several separate three dimensional hydrophone 
arrays sitting on the sea bottom. They are lar ed out in a 
rough grid, each array being approximately 7503 feet from 
the next, so that the sound source being tracked is never 
more than about 5000 feet from the nearest array. The 
arrays are at depths cf roughly 1000 to 1300 feet. 

Each array (see figure 1.1) has four indepei dent hydro- 
phones defining an orthogonal coordinate system. The phones 
are referred to as the x,y, z, and c hydrophones, and are on 

four adjacent vertices of a cube (see figure 1.2'= with sides 
of length D (usually 30 feet). The arrays are linked to 
shore based computers by electronic cable. Th 2 origin of 
the local coordinate system of each array is at the acoustic 
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Figure 1.1 Acoustic Signal Detection by 
Three Dimensional Hydrophonic Array. 
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Figure 1.2 Geometry of Hydrophonic Arrays. 

The sound source is equipped with a clock synchronized 
with the shore based computers and emits a signal at speci- 
fied intervals. The times of receipt of those signals at 
the four hydrophones are recorded^ and the corresponding 
travel times are calculated by subtraction of the signal 
emission time. 

D. APPAEENT POSITIOSS 

The first step in estimating the position of a sound 
source is to do so under the assumption that tha sound wave 
travelled its entire trajectory through water which had a 
constant speed of sound. The result, called the apparent 
position, is obviously erroneous, but is then corrected by 
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the ray tracing procedure (a description of wiich follows 
later). The constant velocity value used is usually the 
estimated speed of sound V at the depth of : he sensing 
array. 

The constant velocity assumption implies that the wave- 
front is an expanding sphere. Therefore the apparent posi- 
tion is calculated using simple spherical eguations 
involving squared distance calculations. Specifically, if 
Tx is the travel time recorded for the x-phone, then V«Tx 
is the distance between the apparent position and the 
x-phone. This distance is also equal to the usual geometric 
distance between the two positions 

( X , Y , Z ) (apparent position) 

and ( D ,-D ,-D ) / 2 (x hydrophone position) . 
E-^uating these two squared distances, and equating their 
counterparts for the other three hydrophones, tie equations 
(1.1) are obtained. 



( X - D/2 )^ + ( Y + D/2 )^ + ( z + D /2 
( X + D/2 ) + ( Y - D/2 )^ + ( Z + D/2 

y 

( X + D/2 ) + ( Y + d/2 ) ^ + ( Z - D/2 ) ^ 

z 

( X + D/2 ) ^ + ( Y + D/2 ) ^ + ( Z + D/2 ) ^ 



Assuming that the times Tx, Ty, Tz, and I; are known, 
(1.1) is a system of four equations in three unknowns X, Y, 
and Z. This overdetermined system will, in general, have no 
exact solution. In fact, even if the time values were 
exactly correct, the system would still have no exact solu- 
tion. This is because the equations correspond to the 
straight line ray paths due to the constant velocity assump- 
tion, whereas the time values correspond to tie true ray 
paths which are not straight due to the actual variation of 
velocity of sound along the ray path. This is a subtle, but 
very important point. 
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so as to 



To throw out one of the equations arbitrarily , 
reduce the system to three equations in three unknowns, is 
to throw away information from one of the hydrophones. The 
psuedo solution currently utilized is to subtract the fourth 
equation from each of the first three, yielding a system of 
three equations in three unknowns which allo#s an exact 
solution involving information from all four hydrophones. 
However that solution will not, in general, satisfy any of 
the original four equations, and is only one of many reaon- 
able ways to choose an approximate solution. 

E. INITIAL ELEVATION ANGLE AND EAI TfiANSIT TIME 

Assuming that a solution (Xa,Ya,Za) has beei determined 
for the apparent position, then the initial angle of eleva- 
tion is just the angle of elevation of that solution, given 
by (1.2). 




The objective is to find an apparent positioi (Xa,Ya, Za) 
such that (1.2) computes an angle which approximates the 
physically correct elevation angle as closely as possible. 
The solution (Xa,Ya,Za) and the times Tx, Ty, Tz and Tc can 
then be used to determine an appropriate value for T, which 
is the ’ray transit time’, or time of arrival of the sound 
wave at the acoustic center of the sensing array. The 
currently employed method uses the proportional relationship 
of equation (1.3), where S and Rc are the ranges from the 
apparent position to the acoustic center and to the c hydro- 
phone respectively, as in (1.4). 
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R 

T 



or 



T 



(1.3) 



R 

c 



T 

c 



V, 



T R 
c 



R 

c 




(1 .4) 



F. RAY IRACIHG 

Whichever method is selected to produce tie apparent 
position (Xa,Ya,Za) , it is transformed into the estimate of 
the true sound source position (X, Y/Z) by the procedure of 
ray tracing. When there is velocity layering ii the water, 
the ray path is no longer a straight line. This is treated 
using repeated applications of Snell’s Law [Ref. 2 p.131], 
starting with the layer of water in which the array sits, 
and backtracking upwards through successive velo: ity layers, 
until the estimated ray transit time T is consumed. 

The layering effect is artificially induced by the limi- 
tation that the speed of sound can be estimated at only a 
finite number of depths, the result of which is commonly 
called the water column. For example, at the tracking range 
studied the speed of sound is measured every 25 feet, 
starting at the depth of 12.5 feet. Hence, for example, the 
third layer from 50 to 75 feet deep is assumed to have a 
constant speed of sound equal to that measured at 62.5 feet. 

The first layer processed [Ref. 3 p.4] is the partial 

layer lying between the array and the deepest la^ er boundary 
that is shallower than the array, with thickness Z1 (see 
figure 1.3). 

The incremental slant range in the first Layer is Si 
given by (1.5), where A1 is the initial elevation angle 
estimate. 
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Figure 1.3 Bay Tracing an Apparent Position. 



Si = z^j sin (Ai) 

The incremental travel time in the first Layer is T1 
given by (1 .6) , where VI is the velocity estimated for the 
layer in which the array is situated. 



( 1 . 6 ) 



Ti = 

The incremental horizontal distance travellel by the ray 
in the first layer is Hi ^iven by (1.7). 



cos(A^) 



(1.7) 



To determine the angle of elevation in the next layer, 
Snell’s Law (1.8) is applied, where V 2 is the speed of sound 
estimated for that layer. 



cos (Aj^) 




When (1.8) is solved 
entry into the next layer, 

cosCa^) = 



cos (A^) 




for the cosine of ti e 
(1.9) is obtained. 
cos(A^) 




( 1 . 8 ) 



angle of 



(1.9) 



The procedure of computing the incremental values of 
slant range, time and horizontal distance are now repeated 
for the second layer. The overall procedure is repeated 
upwards through layers 2,...,n , where n is the first layer 
in which the sum of the incremental travel times exceeds the 
total ray transit time, as in (1.10). 



(Ti . Tj . 



+ T 'N > T 



( 1 . 10 ) 



In the last, uppermost layer the values Tn, Sn, Hn, and 
Zn must be adjusted to compensate for overshooting the total 
time T. The values Hi and Zi (i=1,...,n) are then accumu- 
lated to get (1.11). 



H = 




n 




Now the raytraced position estimate is given by (1.12). 
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X 



( 1 . 12 ) 



■ ' \ " >/V?^ 



Y = ( Y H ) 

3i 



H 



a a 



Z = z 

The sensing array is usually not aligned wi: h the coor- 
dinate system of the overall tracking range. Therefore the 
apparent position, which is in terms of the local array 
coordinates, must be changed by a suitable geoms trie trans- 
formation prior to ray tracing so as to account for the 
angle of tilt at which the array sits on the sea bottom. 
After ray tracing, the position estimate must be again 
transformed to account for rotation of the array about its Z 
axis away from a position which is aligned with the range 
coordinate axes. Finally a simple translation must be 
applied to reference the position estimate to the range 
coordinate system origin vice the acoustic cei ter of the 
array. The end result is a position in terms of the overall 
range coordinate system. These transformations are not 
given here because they are used after the estimation of the 
initial angle and time, and hence do not affect the accuracy 
of those estimates. See [Ref. 3] for further dec ails. 



G. DISCOSSIOH 

If the velocity versus depth relationship is smooth and 
estimated accurately, then the ray tracing procedure is 
surprisingly robust with respect to the thickaess of the 
layers. For example, if velocity is linear versus depth, 
and is known exactly, then the exact hydrophone times, ray 
transit time and initial angles can be computed [Ref. 4] for 
any given sound source position. Then the ray tracing 
procedure, with layers as thick as 25 feet and targets as 
far away as 3000 feet, still estimates positions to within- 
inches of each of the true coordinate values. This seems to 
indicate that errors resulting from position estimation are 
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not due to ths approximation by layers# except pjssibly when 
there are radical changes in the velocity pattern within 
single layers such as frequently occur in layers near the 
water surface. Rather such errors apparently ars due mostly 
to inaccuracies either in the estimation of the speed of 
sound profile itself, or in the initial angle and transit 
time estimates. This study will focus on those srrors which 
are involved in the time values observed at the hydrophones, 
and attempt to produce methods of initial angle and transit 
time estimation which reduce the effects of those errors on 
the overall position estimation procedure. 
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II. C DEBEHT LY EMPLOYED METH ODS 

1. BASIC CONCEPTS 

The method currently used for estimation of an initial 
angle and ray path transit time focuses on the four spher- 
ical equations (1.1) given in Chapter I. A> previously 
discussed, these equations have no exact solution because; 

1. they form an overdetermined system of four equations 
in three unknowns; 

2. the time values recorded at the hydrophones may be 
innacurate due to unknown sources of error in the 
observation process; and 

3. even if the time values were exact, the^ correspond 
to a nonconstant velocity profile, and so will not be 
be correct for the constant velocity geometry (spher- 
ical wavefront) assumed by the equations. 

As previously mentioned, the psuedo solution chosen by 
the current method is to subtract the fourth spherical equa- 
tion from each of the first three, and solve the resulting 
system of three equations in three unknowns. This method 
has the beneficial quality that information is retained from 
all four hydrophones, whereas to just drop the fourth equa- 
tion (or any one of the equations) without the initial 
subtraction would cause complete loss of the information 
from the data recorded by one of the hydrophones . However 
it is important to realize that the solution thus obtained 
does not actually satisfy any of the original four 
equations. 

It should be noted that the development of ta is solution 
in [Eef. 3] is done entirely from a geometrical point of 
view, and does not mention the system of foi r spherical 
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equations. The text of [Ref. 3] does not draw i ttention to 
the fact that the solution developed is just one of many 
plausible choices, none of which will satisfy all four 
spherical constraints. Therefore the solutioi chosen is 
treated as though it were the exact solution, only subject 
to errors in the observed hydrophone time values. However, 
even with exactly correct time values, this currently 
employed method will not yield the true elevation angle and 
ray transit time. This is due to the assumption of a 
constant velocity vice the true nonconstant velocity 

profile. This conflict introduces an automatic bias in the 
initial angle and time estimates currently used for ray 
tracing. 



B. CCMPOTATIONS 



To simplify notation, let (X,Y,Z) be the coordinates of 
the apparent position which was formerly denoted (Xa,Ya,Za). 
Then when the current solution is applied, the first step is 
to subtract the fourth equation from each of the other 
three, which produces the equations (2.1). 

{ X - D/2 - ( X + D/2 ) 

<= X (2.1) 

( Y - D/2 - ( Y + D/2 )“ = ) 

c y 

( Z - D/2 - C Z + D/2 ) 

c z 

The solution to these are easily obtained^ as in (2.2) • 



2 2 9 

X = v (T-T )/2D 
c X ^ 

Y = - tJ ) / 2 D 

2 2 0 

Z = V (T -T )/2D 
c z 



( 2 . 2 ) 



Then the initial elevation angle estimate is (2.3). 



A = 




(2.3) 
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The ray path transit time to the acousti; center is 
(2.4), where R and Rc are as defined in Chapter I by (1.4). 

T = (2.4) 

This method -of determining the apparent position shall 
hereafter be referred to as the ’Navy Method’, or ’NAVY’ for 
short . 



C. ADJUSTMENTS TO THE ORIGINAL SOLUTION 

Experience has shown that the NAVY method produces an 
apparent position estimate which usually can be improved by 
an adjustment which is described in this section,. 

The cosine of the angle between the i-th axis and the 
straight line from the origin out to the apparent position 
is called the i-th direction cosine Ci. It is a fact of 
geometry that the sum of the squares of the three direction 
cosines must equal unity. Therefore the method is now 

adjusted to reflect that constraint. 

The direction cosines used are the, angles nade by the 
ray path at the c hydrophone. Therefore the (X,Y,Z) coordi- 
nates calculated by the original method are first translated 
to coordinates referenced to the c-phone as the temporary 
origin, as in (2.5). 



X = X + D/2 y = y + D/2 Z = Z + D/2 ,o c, 

c c c (2.5) 

Therefore the three direction cosines ars given by 

( 2 . 6 ) . 



c 

X 



X., / V T C 

c c y 



y / V T 
c c 



C 

z 



Z / V T 
c c 



( 2 . 6 ) 



The denominators in (2.6) are all V*Tc beci use that is 
the range from the apparent position to the c hydrophone, as 
estimated by the time from the c hydrophone. Ideally these 



23 



cosines should add to unity when squared. Thera fore if DCC 
is the ’direction cosines correction' factor defined by 
(2.7) , then DCC should be close to one. 



DCC 




(2.7) 



Deviation of DCC from unity is interpreted as an indica- 
tion of receiver timing errors, array malformation or 
invalid data at one or more of the hydrophobes [Ref. 3 
p.C-3], Currently if DCC lies outside tie interval 
(0.98,1.02), the data is discarded as being excessively full 
of error. The direction cosines of the remaining acceptable 
data points are rescaled using (2.8) to assure satisfaction 
of the direction cosines constraint. 



c; = / DCC c- = Cy / DCC c; = C^ / DCC (2.8) 

A corrected set of new coordinates are computed by 
(2.9), still being referenced to the c hydrophone. 



X = VTC Y = VTC Z=VTC 

C CX C cy c cz 



(2.9) 



These are then translated by (2.10) 
referenced to the acoustic center. 



to coordinates 



X = X - D/2 
c 



Y = Y - D/2 
c 



Z = Z - D/2 
c 



( 2 . 10 ) 



This adjusted method of determining the apparent posi- 
tion shall hereafter be referred to as the 'Navy Adjusted 
Method', or ’NAVY_A' for short. 



D. DISCOSSIOH 

When a sound source is within the detectibn range of 
more than one sensing array, each array produces timing 
data. The data from each array can be processed to produce 
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a position estimate. Ideally these multiple estimates of 
position will be in reasonably close agreement. However 
experience with actual tracking data has shown that this is 
not the case in many of the multiple detection opportuni- 
ties. This tendency toward disagreement between multiple 
estimates of the same position is commonly called the cross- 
over, or crosstalk, problem. This problem often occurs when 
the sound source is moving away from the tracking domain of 
one array into the tracking domain of another. This study 
focuses on improvement of the initial angle and time esti- 
mates, which hopefully will help alleviate the crossover 
problem. 

The current choice of a ’best' compromise solution 
appears to be based on reasons of simplifying geometry and 
calculations. These are worthwhile objectives, but do not 
in themselves reflect the need to estimate accurately the 
initial elevation angle and ray transit time- Since there 
exist physically correct values for both the aa gle and the 
time, those values will produce the exact position after ray 
tracing, provided that the velocity -profile is known 

exactly. The desire then is to estimate these true values 

as accurately as possible. 

The question at this point is whether or not the direc- 
tion cosines adjustment causes the estimated apparent posi- 
tion to be closer to the true apparent position. Experience 
has indicated that it does [Bef. 3 p.C-7], However the 

effect of the adjustment can be interpreted in terms of the 
original four spherical equations (1.1). The rescaling of 
the direction cosines given by (2-6), so as to assure that 
their squares add to unity, is equivalent to rescaling the 
quantities in (2.5) so as to assure that their squares add 
to (V«Tc)2. That is exactly the constraint stated by the 
fourth spherical equation of (1.1). So the effect of the 
adjustment is to require that the fourti equation. 
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concerning,' the data at hydr 
This requirement will, in 
three equations are not sa 
that the adjustment often i 
to imply that the fourth e 
than the other three. Or 
faction of one of the eguat 
good compromise solution. 

To summarize, the N 
apparent position suitable 
the direction cosines ad jus 
for reasons not understood 
that position as indicated 
are supported by the result 
However, as will be demons 
in the next section, the DC 
method has an effect which 
the smoothing of timing err 



ophone c, is exactly satisfied. 

general, assure that the other 
tisfied. Since experience shows 
mproves the solution, this seems 
guation is somehow mo: e important 
it may just be that exact satis- 
ions usually assures a reasonably 

AVY method provides a useable 
as input for ray tricing. But 
tment used in the NAVY_A method, 
at this time, appears to improve 
by test results. Ti ose results 
s of this thesis (see Chapter V) . 
trated by the example considered 
C correction factor o: the NAVY_A 

must be something more than just 
ors . 



E. A COHPOTATIONAL EXAMPLE 



For the purposes of illustration and comparison, suppose 
that a 30 foot sensing array is at a depth of 1300 feet, 
that the coordinate system origin is at the array acoustic 
center, and that the array arms are parallel to the coordi- 
nate axes. This implies that the four hydrophones are in 
the positions; 

X : ( 15 , -15 , -15) 

y : (-15 , 15 , -15) 

z : (-15 , -15 , 15) 

c : (-15 , -15 , -15) 

If there is a sound source known to be in position 
( 1000 , 3000 , 900 ) 

then the depth of that source is 1300 - 900 = 400 feet. 
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DEPTH (feet) 



1 




Figure 2. 1 Sample Depth Versus Velocity Profile. 
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Figure 2,1 shows the estimated sound velocity profile 
which was estimated for the data used in the course of this 
thesis. As can be seen, the profile is primarily linear at 
depths greater than 100 feet. The profile at depths greater 
than 100 feet is reasonably approximated by the linear 
relat ionship 

V = 4840.7 + 0.03314 • DEPTH . 

Therefore suppose that in this example problem the 
velocity profile is known exactly, and is given Dy the above 
linear relationship. 

Under these circumstances, witn known linear velocity 
profile, and known sound source location, the exact times of 



arrival of the 


sound wave at 


the f 


our 


hydrophones can be 


computed using 


the methods set 


forth 


in 


Appendix A. Those 


exact times (in 


seconds) are: 








Tx 


; 0.6779686893 




Ty 


; 0.6742257788 


Tz 


; 0.6782243197 




Tc 


: 0.6798324156 . 



The corresponding exact values for the initial elevation 
angle, ray transit time and resulting apparent position are 
also directly computable. Those computations will hereafter 
be be known as the EXACT method, and will produce the 
correct true position after ray tra-cing. The results of the 
EXACT method are given in Table I, along with the corre- 
sponding apparent position estimates produced waen the two 
methods, NAVY and NAVY_A, are applied to the (errorless) 
time values. 

At first glance the differences in Table I might seem 
rather small. However it is important to recall that these 
are produced under the ideal conditions of a very smooth and 
exactly known velocity profile. These idealizations are far 
from the realities of a nonlinear velocity profile which is 
estimated by a procedure involving errors which are unknown 
and probably significant. Such realities mighc well cause 
the differences in Table I to be significantly larger. The 
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TABLE 


! I 










Single Example 
NAVY and NAVY 


Comp 
_A a 


arisen of 
ethods 






Method 


Transit Time 
T Tsecs) 




Elev. Angle 
— 






NAVY 


0. 67526969 




15. 26459 






NAVY_A 


0. 67527027 




15. 26461 






EXACT 


0. 67527043 




15. 27002 








Apparent Position Estimate 






NAVY 


( 1005. 957 , 


3017 


. 874 , 868. 143 


) 




NAVY_A 


( 1006.087 , 


30 18 


.259 , 868.255 


) 




EXACT 


( 1005.966 , 


3017 


.899 , 868.555 


) 


, 



nature and size of those differences remain difficult to 
determine until more is known about the the velocity profile 
estimation errors and their effect on the position esti- 
mating process. In any event even the small differences in 
Table I might be magnified during the ray tracing process 
under the conditions cf a realistic velocity profile. 

The differences illustrate the very important point that 
th e dir ect ion co^nes adjustment causes ch an ges in the esti - 
mate s even whe n the time da_^ is free of all er^or. Hence 
the deviation from 1.0 of the correction factor DCC is not 
just due to array malfunction, receiver timiig error or 
other sources of non-valid data as previously assumed. 

In the example above, the NA7Y_A method produced a 
slightly better time and angle value than the 1 AVY method. 
However, in this same example the NAVY_A method produced an 
apparent position estimate which is slightly farther away 
from the 2XACT answer tnan the position estimated by the 
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NAVY method. This is only one example, and its results 
should not be generalized. However it illustratj s the point 
that apparently the true effect of the adjustment may not be 
well understood. 

Furthermore, the adjustment seems to place heavy 
emphasis on the time value recorded at hydrophone c. 
Therefore the effect of the adjustment may well depend 
largely on the accuracy of that one particular data value, 
which is a relatively unbalanced dependence in the presence 
of data errors. 
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III. PLA NAR BA7EFR0NT MODELS 



A. THE PLANAR WAVEFRONT ASSUMPTION 

When a sound wave travels in constant veloci: y water, it 
expands in the shape of a sphere. If the velocity profile 
is variable instead, but is reasonably well behaved versus 
depth, then the expanding wave is a smooth distortion of a 
spherical surface. In either case, if the wave has 
travelled a long distance when it arrives at a hydrophonic 
array, then that small piece of the wavefront # hich passes 
through the 30 foot cube spanned by the array may be approx- 
imated reasonably by a flat planar surface. This approxima- 
tion is the basis of the planar wavefront models developed 
in this chapter. 

B. PLANE EQUATIONS 

A plane in space is fully defined by idei tif ying any 
point (X0,Y0,Z0) on the plane,, and also a vector (C1,C2,C3) 
of unit length which is perpendicular to that plane. The 
vector is called the unit normal vector for that plane. 
Then any point {X, Y,Z) on that plane must satisfy the equa- 
tion of the plane, namely (3.1). 

Cf ( X - Xq ) + ( Y - Yq ) -H C 3 ( Z - Zq ) = ° ( 3 . 1 J 

The perpendicular distance from the plane t) any point 
(X1,Y1,Z1) not on the plane is the absolute value of (3.2). 

Cl ( Xi - ) + Cj ( Yi - ) * C 3 ( 2^ - Z„ ) 
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C. BASIC flODEL EQOAIIONS 



Consider a coordinate system whose origin is at hydro- 
phone c, and whose axes are aligned with the three array 
arms. Let Cl, C2 and C3 be the direction cosines on the X, 

y and Z axes respectively for the vector from tiie origin to 

the apparent sound source position. Then (C1,C2,C3) is 
itself a vector, of unit length, which is perpendicular to 
the planar wavefront emanating from the soi nd source. 
Therefore (Cl, C2, C3) can be used as the normal vector for 
the wavefront plane. 

In the coordinate system referenced to the c-phone as 
the origin, the acoustic center has coordinates (D,D,D)/2, 
where D is the length of an array arm. When the soundwave 

plane arrives at the acoustic center, it will have the 

equation (3.3) . 



The x-phone has coordinates (D,0,0), and : he distance 
between it and the soundwave plane at the acoustic center is 
(3.4), which then simplifies to (3.5). 



This distance is measured in a direction perpendicular 
to the wavefront plane, and so is measured in the direction 



X + Y + z 




( D/2 - D ) + (D/2-0 ) + (D/2-0 ) 



D ( -Cl + C 2 + C 3 ) / 2 



(3.5) 



of travel of the soundwave. Therefore the same distance is 
also equal to (3.6) . 



V { T - T ) 

X 



(3.6) 
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In (3.6) V is the velocity of sound at the array, Tx is 
the time of arrival of the sound wave at the x-phone, and I 
is the time of arrival of the sound wave at the acoustic 
center. The term ’’distance” is used loosely in (3.5) and 
(3.6) , because these quantities may be negative. The true 
distances are the absolute values of these quantities. 
Since the next step is to equate these two distal ces, it is 
only necessary to show that these two quantities always have 
the same sign. There are two cases to consider, depending 
on whether the first component of the apparent position is 
positive (X>0) , or negative (X<0) . Let (X',0,0) be the 
intersection of the X axis with the wavefront plane as it 
passes through the acoustic center. Then (3.3) can be used 
to solve for X’, namely 

X’ = D (Cl + C2 + C3) / 2 Cl. 

Now consider the case where X>0. Then C1>0 i Iso, and if 
(3.6) is positive, then 



Tx > T 

which implies that the wave arrives at the acoustic center 
before it arrives at the the x hydrophone. Therefore, since 
X>0, the plane at the acoustic center will intersect the X 
axis at a pcint beyond the x hydrophone, or X’>D, and hence 
X’ > D => (Cl + C2 + C3)/ 2 Cl >1 

= > Cl + C2 + C3 >2 Cl (sinc 2 COO) 

= > -C 1 + C2 + C3 > 0 

= > (3.5) is positive . 

A parallel argument applies for the case of X>0, thus 
establishing that (3.5) and (3.6) always have the same sign. 
Equation (3.7) is the result of equating these two 
quantities. 



+ C 2 + C 2 + (2VT/D) - (2VT^/D) 



(3.7) 
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For convenience, let K = 2V/D, and then apply the same 
logic to the distances of the y, z, and c hydrophones from 
the wavefront plane as it passes through the acoustic 
center, to obtain the equations of (3.8). 













- K 


T 




K 


T 

X 


= 0 


-^1 


+ 




- 


s 


- K 


T 


+ 


K 


T 

y 


= 0 


1 

n 


- 




+ 




- K 


T 


+ 


K 


T 

z 


= 0 




+ 




+ 




+ K 


T 


- 


K 


T 

c 


= 0 



The system (3.8) is four equations in four ui knowns, and 
is the planar model's version of the equations (1.1). The 
unknowns in (3.8) are Cl, C2, C3 and T. However there is 
the additional constraint that Cl, C2, and C3 ace direction 
cosines, and therefore the direction cosines constraint 
(3.9) is a fifth equation, creating a a system of five equa- 
tions in four unknowns. 



* 



= 1 



(3.9) 



Generally there is no set of values (C1,C2,:3,T) which 
will satisfy all five equations at once. This is because of 
the realities of a nonplanar wavefront and the presence of 
timing errors. The next section developes a method that 
produces set of values for the unknowns which is intended to 
satisfy those equations reasonably well. 



D. BIHIIJIZATION OF SOM OF SQUARED ERRORS 

Let Ei, (i= 1,2, 3,4) be the value of the left hand side 
of the i-th equation of (3. 8) . Then Ei measures the amount 
of error in the i-th equation caused by the chosen solution. 
It is impossible to have Ei=0 for all i=1,2,3,4. However 
some compromise may be made. Specifically th= compromise 
chosen here is the classic minimization of (3.10), the sun 
of squared errors. 
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(3. 10) 



E 



2 

1 




+ 



E 



2 

3 



minimization is to be done subject to th 2 direction 
cosine constraint (3.9). Application of the Lagrange multi- 
plier technigue [Ref. 5 p. 5 5 ] calls for the minimization of 
(3.11) over all possible choices of Cl, C2, C3, T and 

lambda. 



L 





(3. 11) 



Taking the partial derivative of L with rs spect to T 
yields (3.12) . 

^ T ^ 

= y (-2K E.) = -2K ( E, + E + E + E ) 

Bt i=l ^ 1 2 3 4 (3^^2) 



2 K 



[- 



2K T + K ( T, + T_ + T. 



Equating (3. 12) tc zero and solving for T yij Ids immedi- 
ately the appropriate estimate (3.13) of T, the ray transit 
time. 



— ( Tj + Tj - ) 



(3. 13) 



The partial derivative of L with respect to Cl is 
(3.14). 

3 L 

dc, * ^ ^ 



2 ( El - E2 - E3 + E^ ), - 2 X 

4Cj^ + 2KT + K ( ~ \ 



If (3.13) is used for T in (3.14), then (3.15) results. 




[ ( 4 - 



2 K { T, - T. ) 



] (3. 15) 
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If the same procedure is used for the par:ial deriva- 
tives of L with respect to each of Cl, C2 and C3, and all 
are equated to zero, then the results are (3. 16) . 



( 4 - X ) C . 



2 K ( - T. ) 

4 1 



i = 1,2,3 



(3. 16) 



In order to solve for the Lagrange multiplier lambda, 
square both sides of the three equations of (3.16), and add 
• the resulting equations together. Then use the sum of 

squares constraint (3.9) and solve for Iambi a to yield 
(3.17). 



X 



4 + 



2 K 





(3. 17) 



Substitute (3.17) in (3.16) and simplify to obtain the 
appropriate estimate of Ci, namely (3.18). 



c. 

1 




i = 1,2,3 



(3. 18) 



The choice of sign in (3.18) is determined by the fact 
that Ci is positive if and only if the sound wavs arrives at 
the i-th phone before it arrives at the c-phone, which in 
turn implies that (T4-Ti)>0. 



E. THE LEAST SQUARES METHOD 

In summary, the first alternative method for determining 
an apparent position starts with estimation >f the ray 
transit time to the acoustic center by (3.19). Then the 
apparent position estimates are computed using (3.20). 

T = ( + T 3 - ) / 2 (3, 19) 
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X 



V T ( T - T, ) 

4 1 

V ^ ' ^4 - ’■j 

3 = 1 ^ 



V T { ) 




(3.20) 



V T ( T - ) 

4 3 




This method shall be hereafter referred to a> the 'Least 
Squares Method', or 'L.S. ' for short. The apparent advan- 
tages of the L.S. method are that: 

1. all four hydrophone times have egual # eight in a 
simple expression for the ray transit time T, rather 
than using an expression so heavily dependent on Tc 
as in the NAVY and N AVY_A methods; 

2. the differences of squared time values # hich appear 
in the solutions of the NAVY methods are avoided in 
the L.S. method, thereby lessening tie tendency 
toward computational roundoff problems; 

3. the direction cosines already add to unity when 
squared, requiring no arbitrary adjustmen: ; and 

4. computation of the initial angle allows cancellation 
of several terms, resulting in the simple expression 
(3.21) . 




F. BIAS IN THE LEAST SQUARES METHOD 

Unfortunately a potentially serious problem exists with 
the L.S. method. That concerns the consequeices of the 
assumption of a planar wavefront. The effect is difficult 
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to derive explicitly because it is diificult :o determine 
the true shape of a wavefront corresponding to a nonconstant 
velocity profile. However if it is assumed that a spherical 
assumption is more accurate the planar assumption, then the 
bias can be estimated roughly, and then subtracted from the 



original L.S. position esti 
problem because, as previ 
spherical equations themsel 
Nevertheless, as a ro 
following procedure is use 
position by the L.S. metho 
line distances from that po 
hydrophones. Divide thos 
sound at the array to obta 
these times to recalculate 
L.S. method again. The di 
recalcualted L.S. position 
would be made by the L.S. 
time values which corresp 
sound wave whose source is 
apparent position. Therefo 
a bias vector which can be 
solution. This bias corre 
set forth in the next secti 



mate. This is still a difficult 
ously discussed, the four basic 
ves have no exact solution, 
ugh estimate of the bias, the 
d. First estimate : he apparent 
d. Then calculate the straight 
sition to each of tha four array 
e distances by the velocity of 
in the corresponding times. Use 
the apparent position using the 
fference between the original and 
s roughly measures the error that 
method when it is applied to the 
ond to a spherically spreading 
in the vicinity of the original 
re this difference cai be used as 
subtracted from the original L.S. 
ction is the basis of the method 
on. 



G. THE LEAST SQUARES CORRECTED METHOD 

The second alternative method for estimating an apparent 
position is as follows: 

1. calculate the apparent position PI by using the L.S. 
method; 

2. calculate the distances from that position to the 
four hydrophones, and convert them to times by 
dividing by the speed of sound at the array; 
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3. use the new times to calculate a new apparent posi- 
tion P2 , using the L.S. method; 

4. calculate the difference vector P2-P1 (see figure 
3.1), and suhtract it from the original position PI 
to obtain the corrected position P 

P = PI - { P2 - Pi ) = 2 PI - P2 

5. finally adjust the ray transit time T calculated for 
the original position -Pi, by using the proportional 
transformation: 

T’ = T * R / R1 

where R is the range to the new position ? , and R1 is 
the range to the original position PI. 





Figure 3. 1 Bias Adjustment for the L.S. Method. 

This method shall hereafter be referred to as the 'Least 
Squares Corrected Method', or L.S.C. for short. The proper- 
ties of the L.S.C. method, like those of the NAf Y_A method, 
are not well understood at this time. It is offered only as 
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an alternative which may combine the beneficial properties 
of the l.S. and NAVY methods, namely that it will; 

1. estimate a ray transit time value equally dependent 
on all four data times, thereby smoothing out exces- 
sive error in any one of the time values; and 

2. reflect the spherical wave assumption, believed to 
provide a more accurate description than the planar 
assumption. 

For targets at a range of 3000 feet, the length of the 
error vector P2-P1 varies from 0 to as much as 10 or 11 feet 
(see Table II). The error vector lengths seem to be depen- 
dent cn both azimuth and elevation angles of the target from 
the array. These patterns indicate a potential for further 
investigation to relate estimation errors to sucx variables. 

H. 23AXIHOH LIKELIHOOD CALCOLAIIONS 

One shortcoming of all methods described thus far is 
that none of them can be used to produce an estimate of the 
underlying error in the time data values. The method set 
forth in this section is a first attempt to estimate that 
noise, and is again based on the assumption of a planar 
wavefront. 

Let Ti be the time recorded from by i-th hydrophone. 
Let Di be the true time which, under absolutely error free 
conditions, would have been recorded at the i-th hydrophone. 
An assumption of Gaussian noise is now made, namely that 

Ti = Ui + Ei 

where the Ei are independent identically distributed normal 
random variables with mean zero and variance 32 . Therefore 
the Ti (i=l,2,3,4) are also normally distributed with the 
same variance, but with means Oi (i=1,2,3,4). 
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TABLE II 

Error Vector Lengths for the L.S.C. Method 
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Now let Cj be the j-th direction cosine (j=l,2,3). Then 
geometrically/ using the assumption of a planar wavefront, 
Cj is given by {3.22), where V is the speed of sound at the 
array, and D is the array dimension (30 feet) . 

C, = V ( 0^ - u. ) / D ,3.22) 

Letting U = 04, .then (3.22) can be solved for each Uj in 
terms of 0 and Cj, as in (3.23). 



u . 

3 



U - D C . /V 
3 



(3. 23) 



The probability density of each time value Ti is given 
by (3.24). 

2 

I 

exp (3.24) 









2 S 



Using (3*23) in (3.24) , and multiplying the four densi- 
ties together^ the likelihood function (3.25) is formed. 



L 



0 




I 

1=1 



(Tf - 



U 



+ DC^/V)' 



(3.25) 



In (3.25) C4 has been used for notational ; onvenience, 
and is defined to be zero. Then the log-likelihood function 
is formed by taking natural logs of (3.25), yielding (3.26). 

1 ^ 2 

= -2 In (2 7T) - 4 ln(S) (T^-U+DC^/V) (3.26) 



Since the values of the Ci are to be direction cosines, 
the usual direction cosines constraint must be added to LI 
to form the Lagrangian function L2 given in (3.27). 
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1 



(3.27) 



L 



2 





Li=i 



The objective is to maximize L2 over tie possible 
choices of Cl, C2, C3, U, S and lambda. Ignoring S for the 
moment, maximization of L2 can be acheived by minimization 
of 1 given in {3 •28), 



L = 



2 S 



2 X! ‘ ^ DC^/V)^ 



1=1 



1=1 



(3. 28) 



'ake partial derivatives of L to get (3.29) i nd 
d L 2D 



dc. 

1 



^ (T. - U + DC ./V) - 2 \ C. 



Bl 



r 4 



►U 



y T, - 4« * ^ j: C, 



1=1 



1=1 



(3.30). 

(3.29) 



(3. 30) 



Eguate (3.30) to zero and solve for (J to obtain (3.31), 
tie maximum likelihood estimate for the time at the c hydro- 
phone. 



u = 




(3. 31) 



Eguate the three equations of (3.29) to zero, and 
multiply each equation by Ci respectively, to ob: ain (3.32). 

X^i = i = 1.2,3 (3.32) 

Add the three equations of (3.32) together, and use the 
direction cosine constraint to obtain (3.33). 



X - ^ 



r 3 



y T.c. - u y c. 
1=1 1=1 



D, 



(3. 33} 



Eguation (3,34) results when (3.29) is equated to zero. 
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c, 

1 



T. 

1 



U 



i = 1,2,3 



(3.34) 



V \ D 

n A “ V 



Substitute { 3 . 33 ) into (3.34) to obtain (3.35)/ the 
rnaximum likelihood estimate of the direction cosines, where 
U is the estimate calculated by (3.31). 



c. 

1 



T. - U 
1 



z 

j'l 



T .C . 
3 3 



3 ' 






(3.35) 



As the reader will perhaps have noticed, tie equations 
of (3.35) define each of the unknowns Ci in terms of all 
three unknowns. Such a structure suggests that (3.35) can 
be used as an iteration function. That is, i: reasonable 

initial values are used for the three unknowns in the right 
hand side of (3.35) then new values are produced. Repeat 
the process using the new values until the answer stabilizes 
within acceptable tolerances. Although convergence to the 
correct solution is not guaranteed, the method has never 
failed for the equations of (3.35). Unlike the L.S. method, 

(3.35) does not have any known closed form solution. 

Returning to the standard deviation S, take the partial 
derivative of (3.27) with respect to S to get (3. 36) . 



3^0 1 ^ , 

■■ 2 . = + — - ) T. - u 

ds 2 ,3 A 1 



1=1 



V i ' 



(3.36) 



Multiply (3.36) by and solve for 3^ to get the 

maximum likelihood estimate of the variance, given by 

(3.37) . 



4 



X(T -u) 

1=1 i=l 



(3.37) 



44 



I. 



THE HAXIHaa LIKEIIHOOD PLANAE METHOD 



In summary, the third method for estimation of an 
ai'parent position is as follows; 

1. let U = T4 initially; 

2. use the L.S. method to calculate the initial values 
for Ci, (i=1,2,3), using (&e321) ; 

3. use U and Ci (i=1,2,3) in the right hand side of 
(3.35) to obtain new estimates for Ci (i=l ,2,3) ; 

4. recalculate 0, using (3.31); 

5. reiterate steps 3 and 4 until tne values Ci (i=1,2,3) 
and U converge within acceptable toierancas; 

6. calculate S2, using (3.37) ; 

7. calculate the estimated apparent position relative to 
the c-phone, using (3.3 8) ; 



X 

c 



V U 



Y 

c 



V u 



z 

c 



V u 



(3. 38) 



8 . 



lastly translate 
time estimate to 
acoustic center, 
defined by (1.4) 



this solution and its o rresponding 
a solution and time relative to the 
using (3.39), where R and Rc are as 
in Chapter I. 



X = X + D/2 
c 



Y = Y + D/2 
c 

T = U R / R 

c 



Z = 2 + D/2 

^ (3.39) 



This method shall be hereafter referred to as the 
'Maximum Likelihood Planar Method', or M.l.P. for short. 
Originally the hopes for this method were rather high, espe- 
cially since it was the first method to producs a variance 
estimate. However subsequent experience with the method 
indicates that it probably suffers significantly from at 
least two factors: 

1. the planar wavefront assumption probably builds in a 
position bias as in the case of the L.S. aethod; and 
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2. the variance estimate is inflated since part of the 
noise being measured is due to the inadequacy of the 
planar assumption. 

J. COMPDTATIONAL EXAfiPLE 

for a quick comparison, the three methods ieveloped in 
this chapter are now applied to the example which was used 
in Chapter II. .The EXACT results calculated previously are 
included in Table III for comparison. 





TABLE 


III 






Single Example Comparison of 
Planar Wavefront Methods 




Method 


Transit Time 
T Tsecs) 


Elev. Anqle 

i~Ti§£ir 




L.S. • 


0.67529319 




15.22798 




L. S. C. 


0. 67527149 




15.26372 




M.L. P. 


0. 67529007 




15. 10437 




EXACT 


0. 67527043 




15.27002 






Apparent 


Position Estimate 




L.S. 


{ 1003. 809 


9 


3019.751 , 866. 1> 6 


) 


L.S.C. 


( 1 006. 089 


9 


3018.266 , 868.235 


) 


M.L. P. 


( 997.722 


9 


3023.685 , 859.355 


) 


EXACT 


( 1005. 966 


9 


3017.899 , 868.555 


) 



As can be seen easily, the M.L.P. method fares rather 
poorly in all regards, even worse than the L.S. method. 
This will be confirmed by the evaluations made in Chapter V. 
Also worthy of note is the apparent tendency o: the L.S.C. 



46 



method to correct the L.S. method back toward the exact 
solution. The evaluations of Chapter V will confirm that the 
L.S.C. method almost always yields a better sd lution than 
the original L.S. solution. 
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IV. A SPH ERI CAL WAVER FRO NT MODEL 



A. TBE SPHERICAL WAVEFRONT ASSOMPTION 



All the models developed in Chapter III ire limited 
primarily by the assumption that the wavefront is planar 
upon arrival at the hydrophone array. In this chapter a 
model is developed under the assumption that the wavefront 
is spherical upon arrival at the array. If the sound 
velocity prefile were constant with depth then the spherical 
model would be exact. This is of course not the case, but 
it is suspected that the wavefront is better modelled as a 
sphere than as a plane because that small piece of the wave- 
front which passes through the 30 foot cube spanned by the 
array is locally spherical. That is because e/ ery part of 
that piece travelled through approximately the same regions 
of water, experiencing the same general raytendiig patterns. 

The spherical assumption is accurate if and only if the 
speed of sound is censtant over the ray path, and conseq- 
uently the original four spherical equations apply once 
again. They were: 



( X - D/2 )^ + ( y + D/2 )^ + ( S + D/2 
( X + D/2 )^+ (Y-D/2 )^+ (2+ D/2 
( X + D/2 )^ + ( Y + D/2 )^ + ( Z - D/2 
( X + D/2 )^ + ( y + D/2 )^ + ( Z + D/2 



n o 
= V‘’T" 

X 



2^2 



(4.1) 



= V 



2 2 

= V T 

z 

2 2 

= V T 

c 



It has been stressed previously that there is no exact 
solution (X,Y, Z) satisfying all four equations (:.1), That 
is because the time values on the right hand side correspond 
to the reality of a variable velocity profile. However, if 
the spherical wavefront assumption is to be accurate, then a 
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constant velocity is the assumed case, and any innacuracies 
in the time values are regarded as due to timing errors 
only. Therefore, under the spherical assumption, if the 
true time values Ui (i=1,2,3,4) were known and substituted 
into (4.1), an exact solution to the overdetermined system 
would be realized. In that case a solution to any three of 
the equations would also be equal to that unique exact solu- 
tion. In particular, the NAVY solution of Chapter II would 
be the true solution. Therefore in terms of the coordinate 
system referenced to the c hydrophone, X would be given by 
(4.2) . 



D 

2 



D 



2 2 
( n, - 



However, X is also given by (4.3), wherj 
direction cosine along the X axis of the vector 
hydrophone to the sound source. 



(4.2) 

Cl is the 
from the c 



The time value of interest at the moment is J 4, the time 
at the c hydrophone. Therefore let U = U4 for clarity of 
presentation, and then equate the expressions in (4.2) and 
(4.3) in order to solve for U. If this same Ijgic is also 
applied to the similar expressions for the distances to the 
y and z hydrophones, the results are (4.4). 

= •>/ U“ - (2DUC^/V) + (D/V)^ i = 1,2,3 (4.4) 

The expression (4.4) will be useful in the development 
of the model of this chapter. 



B. LEAST SQUARES MODELS 

A direct approach might be to apply the Ij ast squared 
error technique to the spherical equations, in a manner 



paralleling that used on the four planar eguatLons of the 
L.S. model in Chapter III. However the formulae and equa- 
tions that result are exceedingly complex^ involve fourth 
degree powers of the data values Ti, and have thus far 
defied all solution attempts. Therefore this id 2 a was aban- 
doned in favor of the maximum likelihood approach which 
follows. 



C, MAXIHOa LIKELIHOOD COaPOTATIONS 



As in the M.L.P. model, Gaussian noise is assumed for 
the time data values. Therefore 

Ti = Ui + Ei i = 1,2,3 

where the Ei are independent identically distributed normal 
random variables with mean zero and variance . The 
density of each Ti is therefore (4.5) . 



■r 



(T^) 





i = 1.2, 3, 4 



(4.5) 



To form the likelihood function, multiply the four 
densities together, to obtain (4.6). 



• ■ ( 7 ^) 



exo 






2 Z ^ ^‘4-^’ 



L 2 S“ 1 = 1 



(4.6) 



Form the log likelihood function by taking nitural loga- 
rithms of (4.6) to obtain L 1 of (4.7) . 

r 3 



= -2 ln(27T) - 4 ln(S) - 



y (T.-U.)^ + (T -U)“ 
1-*, 1 1 4 

!.= 1 



(4.7) 



In order to maximize LI with respect to Cl, C2, C3 and 
U, it is sufficient to minimize L2 of (4.8), where a substi- 
tution for Ui has been made using (4.4). 
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= (T-U)^ + £ 



1 = 1 



- V 



U“ - (2DUC^/V) + (D/V) 



(4.8) 



Now add the usual direction cosines constraint and form 
the lagrangian function L of (4.9). 

r 3 



'2 -X 



Vc - 1 

A ^ 



(4.9) 



For rotational convenience/ let Ki be given by (4.10). 



2 o 

= U - (2 D U / V) + (D/V) 



i = 1,2,3 (4. 10) 



Take the partial derivative of L with respect to Ci to 
get (4.11) . 






^2 (T. - K.) 

1 

-2 K, 

1 



■2 D U 
V 



•)-^X=i i. 1,2,3 



Simplify (4.11) and eguate to zero to yield (4.12). 



C. 

1 




i = 1,2,3 



(4. 12) 



Multiply the three equation in (4.12) by Ci (i=1,2,3) 
respectively, and add them together. Then use the 
constraint on the sum of squares of the direction cosines to 
obtain (4.13). 



X 




(4. 13) 



Substitute (4.13) into (4.12) to get the maximum likeli- 
hood estimate for the direction cosines, as in (4.14). 
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c. 

1 



T. 

1 



1,2,3 



(4. 14) 




1 = 



Now take the partial derivative of 
as in (4.15). 




L with n spect to U, 



+ 2 (u - T.) (^* "IS) 

4 



Equate (4.15) to zero and solve for U to obtain the 
maximum likelihood estimate of the time value (4. 16) at the 
c hydrophone. 




(4. 16) 



Finally take the partial derivative of Li in 
respect to S. Equate it to zero and solve to 
the maximum likelihood estimate of the variance. 



(4.7) with 
get (4.17) , 



s 



2 






(4. 17) 



D. SOLUTION BI A HODIFICATION OF NEWTON’S METHOD 

The solutions given by equations (4.14), (4.16), and 

(4.17) once again form a set of equations which would seem 
to be solvable by natural iteration as in the M.l .P. method. 
Unfortunately this time the technique fails to converge. A 
mathematical tool is needed which is stronger chan natural 
iteration. What is used is a modified four dimensional 



52 



version of Newton's .Method [Ref. 5 p.47] to sei rch for the 
roots of a set of four equations. 

The objective is to determine the values Cl, C2, C3 and 
U which satisfy (4,14) and (4.16). For further notational 
convenience define the values M and N as in (4.18) and 
(4.19). 

M 



N 




Given those definitions of M and N, then tie equations 
whose roots are desired can be simplified to (4.20). 



c. 

1 



u 



^4 




1 - N 



Now define error functions as in (4.21). 
the amount of error in each of the equations 
set of values for (C 1 ,C2,C3 , 0) . 

i = 1,2,3 



IhJ se evaluate 
(4.20) for any 

(4.21) 



Let G be the four dimensional coli mn vector 
(gl / 9 2, g3, g4) ' . Finally let GP be the matrix of partial 
derivatives of G, as in (4. 22) . 

Newton's Method in four dimensions says that if X(n) is 
a four dimensional vector holding the current approximate 
roots Cl,C2,C3 and U, then will be an improved 

answer, where X(n+1) is given by (4.23) . 
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(4. 22) 







B^i 


B^i 








S^4 


3^2 






B ^2 






BC 3 








B^3 


B'^3 






BC3 




B^4 




B94 


B^4 






BC 3 


B^ 



Of course it is necessary to calculate the derivatives 
held in the matrix GP in order to use (4.23). Ti ose deriva- 
tives are given in Appendix B. 

Unfortunately when multidimensional versions of Newton's 
Method are applied# there is often a tendency for the method 
to converge slowly# or even diverge. This is because it 
tends to overshoot the best answer for each iteration. To 
alleviate this problem# a modification is made to the 
method. At the end of each Newton iteration# prior to 
proceeding with the next iteration# a Golden Section Search 
[Ref. 6 ] is performed to find the best possibls answer in 
the direction of the new iterative solution. Specifically 
the line in four dimensional space from X(n-1) of the 
previous iteration to X(n) of the present i: eration is 
searched for the best answer. The definition o£ the 'best' 
answer is that point along the search line whir h minimizes 
the sum cf squared error functions, namely (4.24| . 



— n+l 



X 

— n 



Z g- ■ (4. 24) 

i=l 
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The current iterative solution X(n) is thsa given the 
value of the minimizing point resulting from the Golden 
Section Search. Then the next iteration of Newc on’ s Method 
is performed, along with another Golden Section Search to 
find the next iterative solution X(n+1). 

E. THE HAXIHOM LIKELIHOOD HETHOD 

In summary, the fourth and final alternative method to 
estimate apparent positions is as follows: 

1 . let U = T4 initially ; 

2. use the L.S. method (Se321) to initially estimate the 
values of Cj (j=1/2,3); 

3. set X(1) = (C1,C2, C3, 0) ; 

4. initialize the Newton Iteration counter : I = 1 ; 

5. calculate the values Ki, M and N in accordance with 
equations (4.10), (4,18) and (4.19); 

6. calculate the error function vector S (X( I) ) , using 
(4.20) and (4.21) ; 

7. calculate the derivative matrix GP{X(I)), using the 
results in Appendix B; 

8. invert GP(X(I)); 

9. calculate the new Newton estimate X (1+1) from (4,23) ; 

10. perform a Gclden Section Search along the line 
hetween X{I) and X(I+1) to find the point which mini- 
mizes (4. 24) ; 

11. let X(I + 1) be equal to the minimizing point found in 
the previous step; 

12. increase the Newton iteration counter : I = 1+1 ; 

13. reiterate steps 5 through 12 until the values 

X(I) = (Cl,C2,C3,a) 

converge within acceptable tolerances; 

14. calcualte the estimate of S2 using (4,17).. 
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This method shall hereafter be referred to as the 
’.Maximum Likelihood Spherical Method’, or M.L.S. for short. 
As will he seen from the evaluations made in Chapter V, the 
H.L.S. method is apparently the only alternative to consis- 
tently rival the performance of the currently used NAVY_A 
method. It has the additional advantage that it estimates 
the error present in the time data values. 

F. A COflPOlATIONAL EXAMPLE 

This latest method, M.L.S., is now applied to the same 
example considered in Chapters II and III. ? or a guick 
comparison. Table IV lists the results of using all the 
methods. The error vector lengths are the distaaces of each 
position estimate from the EXACT apparent position. 

Since this is only one example, this table is not 
presented for the purpose of any broad conclusioas. However 
it is of interest to note that in this example 

1. the H.L.S. method outperforms all others, including 
the NAVY_A method; and 

2. in many ways the original NAVY method outperforms the 
adjusted NAVY_A method. 

G. VARIANCE ESTIMATION 

In the computational example, the time data values used 
were exact since the methods of ippendix A coi Id be used 
with the known linear velocity profile. Therefore the 
appropriate value for variance in the time vala es would be 
zero. For this error free example, the estimates shown in 
Table V were obtained for the standard deviations of timing 
noise, using the two maximum likelihood methods. 

In the case of this one example, the spherL cal assump- 
tion is apparently an improvement over the planar, since the 

H. L.S. estimate of error is only 4% of the M.L.P. . estimate. 
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TABLE IV 

Single Example Comparison of All Methofls 



1 





Transit Time 


Elev. Angle 


L ength 


of 


M etho d 


T Tsecs) 


1 ri§2§r 




Error"; 


Sector 


NAVY 


0. 67526969 


1 5.26459 




0.412 8 




NAVY_A 


0. 67527027 


1 5.26461 




0.4840 




L. S. 


0. 67529319 


15.22798 




3.739 




L . S. C . 


0. 67527149 


15. 26372 




0.522 




H. L. P. 


0. 67529007 


15. 10437 




13.64) 




11. L. S. 


0. 67527C49 


15.26677 




0.4117 




EXACT 


0. 67527043 


15.27002 










Apparent 


Position 




Estimate 




NAVY 


( 1005. 957 


, 3017.874 


/ 


868. 14 3 


) 


NAVY_A 


( 1 006.087 


, 3018.259 


/ 


868.25 5 


) 


L. S. 


( 1 003.809 


, 3019.751 


9 


866. 12 6 


) 


L . S. C. 


( 1 006.089 


, 3018.266 


9 


868.235 


) 


;i.L. ?. 


( 997. 722 


, 3023.685 


9 


859.33 5 


) 


a. L. S. 


( 1006. 195 


, 3018.190 


9 


868.375 


) 


EXACT 


( 1005. 966 


, 3017.899 


0 


868.55 5 


) 



That is regarded as an improvement because the correct 
answer is zero. The higher M. 1. P. estimate is indicative of 
the inflation due to the planar assumption. 

The tracking range which provided data for this study 
records time values to seven decimal places. Tierefore the 
standard deviation estimated by the M.L.S. method is of 
particular interest, since it indicates errors in the 
seventh decimal place even when there is no ere or present. 
Since the data was actually error free in this ei ample, the 
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TABLE V 

Maximum likelihood Error Estimates 
for the Example Problem 



Model M etho d 
Planar M.L.P. 
Spherical M.L.S. 



1st. Std. D ev il ti on 
5.517 E-6 secs. 

2.191 £■“ 7 s e ' s . 



estimate is a measure of the variation induced bf the spher- 
ical wavefront assumption. A broader discussion of error 
estimation will be undertaken in Chapter VI. 



V. E7ALQA TIC NS OF MODELS 



A. GENEEAL 

There are many sources of errors in the overall hydro- 
phonic tracking problem. These include, among others: 

1. errors in estimation of the sound velocitf profile; 

2. inho mog eneity of the velocity profile over time and 
horizontal displacement; and 

3. possible errors in measuring the positiois, and the 
angles of tilt and rotation for the hydrophonic 
arrays. 

This study focuses on those errors which occur during 
the computations preceding the ray tracing procedure. To 
evaluate the performance of the methods developed in 
Chapters II, III and IV, it is necessary to control strictly 
the ray tracing procedure. Only in that way can the differ- 
ences found between methods be attributed to the differences 
between models, and not to any source of error oi tside those 
methods. 

It was originally hoped that the various methods might 
be compared by applying them to real tracking data. However 
it was found that the overall tracking problem a ad too many 
large sources of error to allow the methods to demonstrate 
any differences. Therefore the methods were cob pared under 
a more tightly controlled simulated environment. 

B. SIHDIATION SCENAEIO 

Two different simulations are used to compare the six 
different methods. Both use a basic scenario similar to 
that of the computational example explored in Chapters II, 
III and IV. That example assumed that; 
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velocity versus depth relationship is linear, and 
n exactly ty; 



V = 4840.7 + 0.03314 • DEPTi ; 



acoustic center of each array is at 
feet ; 

hydrophone arrays are ail level, and 
Z arms are parallel to the respectiv 
of the tracking range, 
these circumstances the methods se 
can be used to compute the exact va 
times, ray transit time and elevation 
emanating from a source at any spe 
efore when the methods are applied to 
resulting estimated positions can be 
rue position. 



a depth of 

their X, Y 
e coordinate 

: forth in 

lues for the 
angle for a 
cified loca- 
those exact 
compared to 



C. SIMULATED ERROR FOR TIMING DATA 

The models developed in this study were designed to 
improve position estimation, especially in the presence of 
errors in the timing data. Lacking any better m^del at this 
time for timing errors, the simulated environment includes 
an assumption of Gaussian errors for the hydrophone times. 
Therefore realistic timing data can be simulate! by adding 
to each exact time value a random quantity of normally 
distributed error. The mean of the error is assumed to be 
zero. The variance was estimated from real tracking data, 
using the variance estimating property of the H.L.S. model. 
The data from one tracking run was used, involving six 
hydrophone arrays and 733 position estimates (sea Chapter VI 
for data selection details) . Each position from the 
tracking run produces one estimate of the variance. The 

variance value chosen for use in the simulations was the 
median of the 733 variance estimates produced bp the M.L.S. 
method. That value was 
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S2 = 9. 1204 E-12 secs2 . 

That is the same as a standarad deviation of 

S = 3. 02 E-6 secs . 

Each of the two types of simulations was rui four sepa- 
rate times. Each run was done with a different specified 
variation. Those four distinct error conditions are delin- 

eated in Table VI . 



I 

TABLE VI 

Simulation Error Levels 



RUN 


LEVEL 


VARIANCE 


STD. 


DEVIAT ION 


1 


zero 


0.0 


0.0 




2 


low 


9.12 E-14 


3.02 


E-7 


3 


medium 


9. 12 E-12 


3.02 


E-6 


4 


high 


9.1 2 E-10 


3.02 


E-5 



D. SINGLE AEBAY SIMOIATION 

In the first simulation, the intent was to compare the 
methods pairwise, so as to determine which method is more 
likely to produce the more accurate estimate of a given 
sound source position. One thousand positions were chosen 
at random. Each position was 3000 feet from the array. The 
positions were uniformly spread over the surface of a sphere 
of radius 3000 feet, centered at the array, truncated above 
by the water surface (depth 0) and below by the depth of the 
array (1300 feet). The methods set forth in Appendix C were 
used to assure that the random positions selected were 
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uniformly distributed over the surface area of the truncated 
hemisphere. 

Each of the 1000 randomly chosen source positions was 
then processed as follows: 

1. calculate the exact hydrophone times, using the 
methods of Appendix A; 

2. add to each of the four exact times a random value of 
error at the specified level (zero, low, medium or 
high) ; 

3. apply each of the six methods to the hydrophone 
times, generating six different apparent positions; 

4. apply the ray tracing procedure to each of the six 
positions, using layers that are 25 feet thick, and 
utilizing the known linear velocity profile, thereby 
producing six different estimates of the sound source 
location ; 

5. compare the six different estimates pairwise to see 
which method in each pair produced the estimate 
closest to the true sound source location. 

The comparison being made is that one method is consid- 
ered preferrable to the other if it more frequently produces 
the more accurate estimate. 

The layer thickness of 25 feet was selected order to 
simulate the actual procedure at the tracking range which 
provided data for this study. However, as discussed in 
Chapter I, when the velocity profile is smooth and known 
exactly, the process is very robust with respect to the 
thickness used. The thickness values 1, 10, and 25 were all 
tried, with virtually no changes in any of the comparisons 
between methods. 
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E. SINGLE ABRAY SIMOIATION RESaLTS 



Tables VII and VIII contain the results of the sin^'le 
array simulation. Each tabled entry represents : he fraction 
of time that the method of that row produced a better esti- 
mate than the method of that column. For example, in Table 
VII, the L.S. method outperformed the M.L.S. method in only 
5.8% of the 1000 trials with low error values. 

In a one tailed test that one method is o etter than 
another, these binomial proportions are significant at the 

0.0 5 level if they exceed 0.526. Symmetrically, one method 
is significantly worse than another at the 0.05 level if the 
proportion is less than 0.474. For the 0.01 level the 
corresponding critical values are 0.537 and 0.463 
respectively. 

The results indicate that: 

1. as previously claimed, the MAVY_A method usually 
outperforms the unadjusted NAVY method; of particular 
interest is the case of zero error whir h actually 
compares the relative ability of each method to 
produce the exact answer when given the exact times; 
in those cases the NAVY_A method does extremely well 
against the NAVY method; 

2. under all error conditions the most successful 
performer is the H.L.S. method, since it always has a 
favorable (greater than 0.5) comparison fraction 
against all other methods; the H.L.S. fractions vary 
little over the four error levels; 

3. under all error conditions, the worst performer is 
always the M.L.P. method; 

4. spherical methods consistently outperform planar 
methods; and 

5. increased error levels tend to lessen the distinction 
between methods; in the zero error case comparisons 
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TABLE VII 

Single Array Simulation Results 
for Lower Error Levels 







ERROR 


LEVEL 


: ZERO 








NAVY 


NAVY_A L 


. S. 


L.S.C. 


M.L.P. 


M.L. S. 


NAVY 




. oil 


.950 


.6 76 


:987 


. 420 


MAVY_A 


.98 9 




-950 


-728 


-987 


. 434 


L. S. 


.05 0 


. 050 




.050 


.987 


.057 


L • S« C • 


.32 4 


. 272 


.950 




.987 


. 346 


M. L. P. 


.013 


. 013 


.01 3 


.013 




.013 


H.L.S. 


.580 


. 566 


.943 


.654 


.987 








ERROR 


LEVEL 


: LOW 








NAVY 


NAVY_A L 


.S. 


L.S.C. 


M.L.P. 


H.L.S, 


NAVY 




. 165 


.952 


.648 


.987 


.396 


NAVY_A 


.835 




.952 


.693 


.987 


-405 


L. S. 


.048 


. 048 




.048 


.987 


.058 


L. S. C. 


.352 


. 307 


. 952 




.987 


.369 


Jl.L. P. 


.013 


. 013 


.013 


.013 




.013 


a. I. s. 


.604 


. 595 


. 942 


.631 


. 987 




are 


in the 


interval 


(0.01 


,0.99) , 


while 


in the h 


error case 


that interval 


is narrowed con 


3 iderably 



(0. 37,0 .63) . 

F. DOUBLE ARRAY SIHDIATIOH 

In the second simulation, the intent again was to 
compare the methods pairwise, this time determining which 
method is mere likely to produce positions which agree more 
closely in the two array crossover problem. This is not the 
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TABLE 7III 

Single Array Simulation Results 
for Higner Error Levels 







ERROR 


LEVEL 


: MEDIUM 








NAVY 


NAVY_A 


L.S. 


L, S. C. 


M.L.P. 


M.L.S. 


NAVY 




.464 


.819 


.559 


.987 


. 435 


NAVY_A 


.546 




.821 


.565 


.987 


. 437 


L. S. 


.181 


. 179 




.181 


.971 


. 177 


L • S. C • 


.44 1 


. 434 


.819 




.987 


. 450 


M.L. P. 


.013 


. 013 


.029 


.0 13 




.013 


M.l. S. 


.56 5 


. 563 


.823 


.550 


.987 








ERROR 


LEVEL 


; HIGH 








NAVY 


NAVY_A 


L.S. 


L.S. C. 


M.L.P. 


M.L.S. 


NAVY 




. 486 


.544 


.439 


. 562 


. 431 


N AVY_A 


.514 




.546 


.516 


.559 


.431 


L . S. 


.456 


. 454 




.457 


.557 


.400 


L. S. C. 


.56 1 


. 484 


.543 




.562 


. 427 


M.L. P. 


.438 


. 441 


.‘443 


.438 




.373 


M.L.S. 


.56 9 


. 569 


.600 


.5 73 


.627 




same question as 


that addressed 


by the single array Simula- 


tion. The 


two 


estimates 


produced by any 


one ma 


thod may be 


very close 


to each other. 


and yet be far 


away from the true 



position. 

For the double array scenario two arrays are used, sepa- 
rated by 7500 feet, both at depths of 1300 feet.. The arms 
of each array are parallel to the corresponding coordinate 
system axes. Once again 1000 positions were raniomly gener- 
ated in a uniform manner, this time over a 3 dimensional box 
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Figure 5.1 Double Array Simulation Configuration. 

running crossways between the two arrays. The box is 1300 
feet deep, 5000 feet long and 1000 feet across (see figure 
5.1) . 

Each of the 1000 randomly chosen sound source positions 
in the two array simulation were processed as follows: 

1. calculate the exact hydrophone times for the first 
array using the methods set forth in Appeidix A; 
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2. add to the exact times some random error at the spec- 
ified level; 

3. Repeat steps 1 and 2 for the second array, using a 
new set of random error values; 

4. apply each of the six methods to the time values from 

both arrays, producing six differen: pairs of 

apparent positions; 

5. apply the ray tracing procedure to bo: h apparent 
positions in each of the six pairs, utilizing the 
known linear velocity profile, producing six pairs of 
estimated sound source positions; 

6. for each of the six pairs of positions, cilculate the 
distance between the two positions in the pair; 

7. make pairwise comparisons of the six different 

distances, to see which method in each pair exhibits 
the closest agreement between its two position esti- 
mates. 

This time the comparison being made is that one method 
is considered ‘better than another if the positions it 
produces agree more closely more often than those of the 
other method. 



G. DOUBLE ARRAY SIMULATION RESULTS 

Tables IX and X contain the results of the 1 ouble array 
simulation. Each tabled entry represents the fraction of 

time that the method of that row produced a pair of esti- 
mates which were in closer agreement than the estimates 
produced by the method of that column. Significance 
criteria for these proportions are the same as for the 
single array simulation. 

These results are similar to those of the single array 
simulation, because they indicate that: 
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TABLE IX 

Double Array Simulatiou Results 
for Lower Error Levels 



ERROR LEVEL : ZERO 





NAVY 


NAVY_A L 


.3. 


L. 3. C. 


M.L.P. 


M.L. 3. 


NAVY 




. 366 


.989 


.726 


.968 


.212 


NAVY_A 


.634 




.989 


.734 


.968 


.212 


L. S. 


.01 1 


. on 




.011 


.555 


.014 


L • SmCm 


.27 4 


. 266 


.989 




.968 


. 193 


M.L. P. 


.032 


. 032 


.445 


.0 32 




. 033 


I1.L.S. 


.788 


.788 


.986 


,807 


.967 








ERROR 


LEVEL : 


LOW 








NAVY 


NAVY_A L 


.3. 


L.3 .C . 


M.L.P. 


M.L. 3. 


NAVY 




.448 


.989 


.661 


.952 


.310 


NAVY_A 


.552 




.989 


.684 


.952 


. 312 


L. S.. 


.01 1 


.011 




.011 


.560 


.014 


L. S.C. 


.339 


.316 


.989 




.952 


.254 


M.L. P. 


.04 8 


.048 


.440 


.048 




.038 


M.L.S. 


.690 


. 688 


.986 


.746 


.962 




1 . the 


NAVY, A 


method outperforms the 


original 


unad j usted 


NAVY 


me t h od 


f although the 


difference is not signifi- 


cant 


in t he 


higher error level cases; 




2 . the 


most successful 


performer 


is consis 


tently the 


M.L, 


S. method; 










3. spherical 


me thods 


almost 


always 


outperf 


orm planar 



methods; 

4. increased error levels tend to lessen the distinction 
between methods. 
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TABLE X 



Double Array Simulation Results 
for Higner Error Levels 







ERROR 


LEVEL 


: MEDIUM 








NAVY 


NAVY_A L 


.S. 


L . S . C . 


M.L. P. 


M.L. S 


NAVY 




. 498 


.84 1 


.513 


.808 


. 443 


NAVY_A 


.502 




.842 


.542 


.809 


. 440 


L.S. 


.159 


. 158 




.159 


.525 


. 160 


L. S. C. 


.487 


. 458 


.841 




.808 


. 433 


M.L. P. 


.192 


. 191 


.475 


.192 




. 144 


M.I. S. 


.557 


. 560 


.840 


.567 


. 856 








ERROR 


LEV EL 


: HIGH 








NAVY 


NAVY_A L 


.S. 


L . S . C . 


M.L.P. 


M . L . S 


NAVY 




. 487 


.556 


.535 


.479 


.442 


NAVY_A 


.513 




.55 8 


.499 


.481 


. 440 


L.S. 


.444 


. 442 




.445 


.451 


.422 


L. S. C. 


.46 5 


.501 


.555 




.4 80 


.438 


M.L. P. 


.52 1 


.519 


.549 


.520 




.428 


M.L. S. 


.558 


.560 


.578 


.562 


.572 





There are however some indications from these results 
which are different from those of the single array simula- 
tion, such as: 

1. the worst performer was consistently the L.S. method, 
rather than the M-L. P. method; 

2. in the high error case the H.L.P. method Is equal to, 
or even marginally better than, every other method 
except the M.L.S. method; 



3. the performance of the H.L.S. method i> noticeatly 
better under error free conditions. 

The last two inferences are perhaps the most inter- 
esting. The maximum likelihood approach was intended to 

estimate and account for Saussian errors in the timing data 
values. Hence it is really not very surprisii g that the 
H.L.S. method does well with error prone data. But it is 

interesting to note that the M.L.P. method also does well 
under high error levels, even though it probably suffers 
from a bias due to the planar assumption. 

Of even greater interest is that the H.L.S. method seems 
to be at its best when compared to other methods under error 
free conditions. This result was unexpected, and indicates 
that the H.L.S. method not only handles timing errors well, 
as was intended, but apparently also does an evei better job 
of approximating the elusive exact solution to the original 
four spherical eguations of (4.1) and (2.1) when the exact 
time values are available. 

H. llflITATIONS ON I NTEEPEETATION OF EESOLTS 

The results of both simulations seem to imply that the 
H.L.S. method outperforms the currently used NAVY_A method 
on any randomly chosen sound source position, with or 
without timing errors. These are encouraging results. 
Nevertheless the reader is cautioned that these tests are 
just simulations, and like all simulations, must make 
assumptions which cannot fully reflect the reali: y of actual 
hydrophonic tracking conditions. The most important assump- 
tions made for these simulations are: 

1. the sound velocity profile is known exactly, and is 
linear; 

2. the errors in the timing values from any hydrophone 
are normally distributed with mean zero, and are 
independent of the noise in any other hydrophone; and 
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3 . 



the parameters of the simulation were fixed; for 
example the single array simulation used a fixed 
range of 3000 feet, and both simulations used arrays 
at 1300 feet of depth, with fixed orientations to the 
range coordinate system. 

The first assumption is probably of little oonseguence. 
It not only greatly facilitates computations, but also helps 
to isolate the initial angle and time estimation problem 
from the unrelated errors involved in the velocity profile 
estimaticn procedure. 

The second assumption is somewhat more troublesome. 
Errors may not be Gaussian at all, or if they are, the mean 
may not be zero. Unfortunately each positioi estimation 
involved only four equations, and therefore did not allow 
for estimation of more than four paramenters. Taerefore the 
error mean, being a fifth parameter, could not be estimated. 
Also the errors of any one hydrophone may very well not be 
independent of the errors of the other three pha nes on that 
array. Fortunately these concerns are offset somewhat by 
the results of the double array simulation, wherein it was 
found that the M.L.S. method was at its best when there was 
no noise at all. 

Also the error type and level may depeni on other 
factors, such as the target's range, elevation and azimuth 
angles from the array. This highlights the concerns of the 
third assumption. There is considerable room here for 
future work concerning the dependency of results on such 
complicating factors. 

Lastly it should be pointed out that the simulations 
make comparisons only on the binomial basis of better versus 
worse in 1000 trials. The magnitudes of the actual differ- 
ences are ignored. It is possible, though perhaps unlikely, 
that while one method marginally outperforms a second method 
in most trials, in all the remaining trials the : irst method 
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There is also room here for 



is much 
further 



worse than the second, 
work. 
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VI. COHC IDSI ONS AND 

A. ESIIHATION OF TIMING DATA VAfilATION 



One of the primary purposes of this study ras to esti- 
mate the amount of variability in the timing data being 
recorded during actual tracking runs. This problem was 
addressed by the M.L. P. and M.L.S. maximum likelihood 
models. The M.L.P. model, as previously discussed, suffered 
from a bias due to the planar assumption, was the poorest 
estimator of positions among all the models, and produced an 
inflated variance estimate. Therefore the spherical model 
M.L.S. is used to estimate the data variance. 

The variability that is measured by the M.L.S. model is 
made up of three components. First there is ; he variance 
induced by the spherical assumption. Then there is the 
variance caused by the seven decimal accuracy used when 
recording the data. Finally there are the errors inherent 
in the physical process, due to such factors as hydrophone 
variability or malfunction, local distortions of the sound 
wave, and inexactness of the water column which estimates 
the speed of sound profile. The last two sourr es of error 
together make up the variability that is involved in the 
time values which are ultimately used in position estima- 
tion, and is therefore the variation that is to be 
estimated. 



If it is assumed that the variability indi ced by the 
spherical assumption is independent of the data variability, 
then the M.L.S. variance estimate is the sum of those two 
variances, or 



cr 
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mis 




+ 




me 
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Therefore the data variance can be estimated by first 
estimating the variability induced by the spherical assump- 
tion, and then subtracting it from the M.L.S. estimate of 
the variation. 

The M.L.S. estimate was obtained by applying the M.L.S. 
method to the data from a tracking run at the Nanocse 
torpedo tracking range on May 6, 1980. That run involved 
position estimation by several different hydrophone arrays. 
The run made several thousand position estimates, 733 of 
which were at depths of 100 feet or more and involved 
targets not more than 4700 feet from the sensing array. The 
depth limitation was imposed to avoid the excessive compli- 
cations caused by the radical changes in the velocity 
profile above that depth (see figure 2. 1) . The maximum 
range limitation imitates the data validation procedure at 
the tracking range, where positions farther than 4700 feet 
from the array are discarded. 

The tracking data yielded the following range of esti- 
mates for the standard deviation of the data noise, using 
the M.L.S. model. 



MAXIMUM VALUE 


2-89 


E-5 


secs. 


0.95 QUANTILE 


1. 18 


E-5 


secs. 


MEDIAN VALUE 


3. 02 


E-6 


secs. 


0.05 QUANTILE 


4. 1 1 


E-7 


secs . 


MINIMUM VALUE 


3. 13 


E-8 


secs. 


Eor an overall estimate 
used, so that: 


of the noise. 


the 


^mls 


( 3.02 E-6 


) 2 





The variability induced by the spherical assumption was 
estimated by applying the M.L.S. procedure to perfectly 
noiseless data in the idealized environment of the single 
array simulation of Chapter V. This was done for targets at 
ranges of 1500 to 450C feet, at 500 feet increments, with 
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1000 randomely chosen targets at each range 
are collected in Table XI . 



The results 



TABLE XI 

Inflation of Error Estimates Induced 
by the Spherical Model 

Standard Deviation Estimates from Exact Timj Values 
with Known Linear Velocity Profile 



RANGE 


MINIMUM 


Q (.05) 


MEDIAN 


Q (.95) 


MA( IMUM 


1500 


2.72E- 10 


2.84E-8 


2. 08E-7 


3. 13E-7 


3. 29E-7 


2000 


2.12S- 9 


5.88E-8 


2. 26E-7 


3. 14E-7 


3.3 OE-7 


2500 


3.90E-8 


1 .09E-7 


2. 34E-7 


3. 14E-7 


3.3 IE-7 


3000 


8.35E-8 


1 .50 E-7 


2.37E-7 


3.15E-7 


3.2 4E-7 


3500 


1.21E-7 


1.59E-7 


2.37E-7 


3.13E-7 


3. 24 E-7 


4000 


1 .47E- 7 


1 .65E-7 


2.39E-7 


3. 16E-7 


3.2 5E-7 


4500 


1 .46 E-7 


1 .69E-7 


2. 42E-7 


3.14E-7 


3. 25E-7 



Table XI shows that the median and maximun inflation 
values are reasonably independent of target range. The 
minimum values vary somewhat/ but only for range values 
below 3000 feet. This represents a very stable situation 
overall. Therefore the inflation due to the spherical 
assumption is estimated by the median value at a range of 
3000 feet/ namely: 



(J 



sph 



( 2.37 E-7 ) 



Combining these two estimates, the variance estimated 
for the timing data is 
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2 




CJ 



ml s 



( 3.02 E-6 ) 
( 3.01 1 E-6 



^ sph 

- ( 2.37 S-7 ) 2 



As can be seen/ the error induced by the spherical model 
is less than 10?? of the H.L.S. error estimate. Therefore 
when it is accounted for by subtraction from the H.L.S. 
variation estimate, the final variance estimite changes 
little. 

The estimated value indicates a standard error in the 
6th decimal place. With a typical speed of sound value of 
4880 feet per second, this represents a position differen- 
tial of about 

4880 • 3.011E-6 = 0.015 feet . 

This estimate is guite low, indicating that the time values 
being recorded are sufficiently accurate. 

There is considerable opportunity for addi: ional work 
determining the relationship, if any, between the time vari- 
ance and other factors such as angles of elevation and 
azimuth of the target from the array. 

There is also the problem that the time / ariation is 
likely to be array dependent. For example, consider figure 
6.1 , wherein the standard error estimates from the actual 
tracking run are plotted versus the range of the target. 
The plot does not indicate that there is any simple rela- 
tionship between range and error level. However the plot 
does show a bunching pattern. When the error estimates are 
plotted separately for each array, then the bunching pattern 
becomes clearly associated with the individual arrays. 
Consider figures 6.2 and 6.3, where the separate plots have 
been made for four different arrays. It is still not clear 
from these plots whether the principle effect is due to the 
individual arrays, or the ranges of the targets . However 
some level of array dependency seems likely, indicating a 
need for additional investigation. 
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DATA FROM ALL ARRAYS 





Figure 6.1 Error Estimation Versus Range of Target. 

B. CHOICE OF METHOD 

Clearly all indications are that the plani r wavefront 
models, L.S. and M.L.P. are not candidates for use as posi- 
tion estimators. Furthermore the hybrid model L.S.C. is an 
interesting improvement, but never really performs well 
enough ccmpared to the M.L. S. and NAVY models. 

The original NAVY model is usually outperformed by the 
adjusted NAVY_A method. However the differences are not 
always significant. 

The spherical model M.L.S., on the other haid, consis- 
tently outperformed all other methods during the simulated 
evaluations. It wculd seem that M.L.S. is tie model of 
choice. It does the best job of handling normally distrib- 
uted errors in the data. But that is not tie strongest 
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1 



DATA FROM ARRAY NO. 7 




DATA FROM ARRAY NO. 14 




Figure 6.2 Error Estimation for Arrays 7 and 14. 
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DATA FROM ARRAY NO. 56 




DATA FROU ARRAY NO. 57 







Figure 6.3 Error Estimation for Arrays 56 and 57. 
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use. 



and surprising. 



argument for its use. A 
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times are used, the M.L.S. 
usually produce the most 
apparent position. This i 
that the sensitive ray tra 
minimally ty apparent posit 
There are nevertheless 
should be considered befo 
wholeheartedly. The first 
namely that these conclus 
idealized conditions of t 
second caution, also previo 
magnitudes of the differenc 
teen ignored. It is con 
always produces a better es 
ence between any two positi 
Lastly there is the 
involves a complicated i 
considerable computer time 
procedure for use with ' 
execution of tracking runs 
NAVY_A method currently in 
to its simple computations. 

However, for post run 
calibration of the hydroph 
recommended as being a raor 
estimator than those method 



more important, 
hat when the exact, error free 
apparent position e; timate will 
accurate estimate of the true 
s the desired overall result, so 
cing procedure will be affected 
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several notes of caution which 
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ions were arrived at under the 
he simulations scenarios. The 
usly stressed, is tha: the actual 
es between position estimates has 
ceivable that while one method 
timate than another, the differ- 
on estimates is acceptably small, 
caution that the M.L.S. model 
terative procedure which uses 
. It is probably too slow a 
real time' analysis during the 
. For real time tracking the 
use is probably preferrable due 



analysis , 
one arrays, 
e exact and 
s currently 



and also .oossibly for 
the M.L.S method is 
more robi st position 
in use. 



C. EECOMMENDATIONS FOR FOTDRE IHQCJIRY 

Several recommendations have already been mi de for work 
needed to estimate the effect of suitable independent vari- 
ables on both timing errors and the bias in certain methods. 
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In addition there exist at least two other areas for 
possibly fruitful investigation. 

The first area concerns the interplay betw 2 en methods. 
Specifically^ the binomial comparisons of Chapter V show 
that even the worst methods are better than each of the 
other methods at least part of the time. Hence it is 

possible that the best method overall would be a suitable 
combination of methods, wherein each method is used where it 
is most effective. For example, even though the M.L.S. 
method has been indicated as the best method for any 
randomly selected position, it may be consisteatly outper- 
formed by another method under certain circumstances, such 
as extremely high or low elevation angles. 

The second area for possible work addresses the question 
of how to next improve upon the existing models. It is 
herein suggested that the next improvement in modelling 
would be a method which is based on a linear velocity 
assumption. As figure 2.1 demonstrates, a linear velocity 
profile is a reasonable approximation for most depths. This 
would be the next logical step above the constant velocity 
assumption which is associated with the spherical models. 
Most of the mathematical basis for such a model in contained 
in Appendix A. Possibly a suitable set of equations could 
be developed involving the hydrophone times and reflecting 
the linear velocity assumption. If so, the least squares or 
maximum likelihood techniques might provide useful results. 
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APPENDIX A 

LIHEAB VELOCITY PROFILE THEORY 



All computations in this study are cade under the 
assumption that the sound velocity is directly related to 
depth ’in a linear manner, and is known exactly. Under these 
circumstances many closed form results, not otherwise avail- 
able, can he obtained and used in those computations. 

Suppose that the velocity profile is given bf 

V = VO + V 1 * z 

where VO and VI are known constants, and z is the depth 
variable, measured down from the water's surface. In this 
case it is known [Hef. 4] that the path of a soi nd ray from 
a signal source to a hydrophone will have the shape of an 
arc of a circle. The center of that circle will be some- 

where above the surface of the water. The vertical place- 
ment of that center is determined by tne value of z at which 
the speed of sound eguals zero (see figure A.1). Although 
that depth is negative and is not really a depth at all, it 
nevertheless has geometric meaning. 

Consider the vertical plane containing the circle 
center, the sound source and the hydrophone. Let h be the 
variable which measures the horizontal position in that 
plane. Let (h, z )= (a 1,a2) be the position of the hydrophone, 
and let (h, z) = (p 1, p2) be the position of the sound source. 
Then C2 given by {A. 1) is the Z coordinate of the circle 
center. 




What must be found is 21, the h coordinate of the circle 
center, and r, the radius of the circle. To solve for these 
values, the equation of the circle is used, evaluated at the 
the two known locations, yielding (A. 2) . 
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Figure A.1 Circular Bay Paths for a Linear Velocity Profile. 



( - Cj ) 



( p, - 



( ) 



= r 



(A. 2) 



( P. 



) 



The left hand sides of the two eguations of (A. 2) can be 
et^uated, and solved for Cl, leading to (A. 3). 

(Pj + a^) (p, - a^) 



2 (p - a, ) 2 2 2 



(A. 3) 



When this value of Cl is substituted intj the first 
equation of (A. 2), then the result is given by (A. 4). 
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The circular arclength between the hydroph) ne and the 
sound source is easily computed. Unfortunately the velocity 
of sound does not stay constant along that arc. Therefore 
in order to determine the amount of time (T) r eguired for 
the sound ray to travel the ray path, the effect of the 
velocity must be integrated along the arc. This is dene by 
(A. 5), where S is the arclength between the two points, and 
V (s) is the speed of sound as a function of position on the 
arc . 

(A. 5) 




In (A. 5) 
arclength of s 
to arclength s 
be equivalent t 



the sound source position corresponds to an 
= 0, and the hydrophone position correponds 
= S. In [Ref. 4] this integral is shown to 
o (A . 6) . 



T 




{A. 6) 



In (A. 6) AO is the angle of elevation of the ray path at 
the sound source, and Al is the elevation ai gle at the 
hydrophone. The antiderivative in (A. 7) can be used to 
solve {A. 6) , leading to the ray path transit time expression 
in (A. 8) . 



/■ 



da 

cos ( a ) 



In 




+ s in ( a ) 
cos ( a ) 



(A. 7) 



T 




■ /cos ( )' 



'1 + sin(A )\ 



/ + sin (Aj ) 



■; 



(A. 8) 



84 



If the elevation angle at any point along the arc is 
denoted hy A, then (A. 9) relates the angle to thi derivative 
of z with respect to h along the ray path. 



tan (A) 



dz 



(A.S) 



Implicit differentiation of the equation of the circle 
yields (A. 10) as another expression for the same derivative. 

h - C, 



dz 

dh 



(A. 10) 



Equating these two derivatives leads to (1.11) which 
relates the elevation angle and the position (h,z) on the 
arc. 

h - C. 



tan (A) 



1 



(A. 11) 



From (A. 11) a simple geometric argument produces the 
equations in (A. 12) . 



z - c. 



cos ( A ) 



sin ( A ) 



(A. 12) 



First let A, h and z be equal to (Al,a1,a2) in (A. 12), 
and then let them equal (A0,p1,p2). Then substitute both 
expressions into (A.S). The result is (A. 13), the desired 
expression for the ray path transit time in t 2 rms of the 
positions of the sound source and the hydrophone. 



T 




/^2 Pi ■ 

V2 - ^2j\^ ^ 



(A. 13) 



To summarize, if a ray path is to terminate i t the three 
dimensional position (X1,Y1,Z1), and the sound source is at 
(X2,Y2,Z2), then perform the following steps in order to 
calculate the exact ray path time and elevatioi angle that 
would correspond to a linear velocity profile: 
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1. use (A. 14) to translate the three dimensional posi- 
tions to the two dimensional coordinates of the 
previously described vertical plane, with the origin 
at the water surface directly above the end of the 
ray path; 



= '\J - Y^) ^ (A. 14) 



2. calculate C2, Cl and r using (A.1), (A. 3), and (A. 4); 

3. calculate the exact ray path transit time T using 
(A. 13); and 

4. calculate the exact elevation angle A at the hydro- 
phone, using (A. 15). 




A 



arccos 




r 




(A. 15) 
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appendix b 

PAETIAl DEEI7AIIVE FOEMULAE FOE NESTON’ S METHOE 

The central tool used hy Newton's .'lethod in the develop- 
ment of the M.L.S. model in Chapter IV is the matrix GP. It 
is the matrix of first order derivatives of the error 
expressions g2, g3 and g4. Those derivat. ves are set 

forth in this appendix. 

If X = (C1^C2,C3,U) , then the error functions are 

defined by (B. 1) , 

T . - /k“ 

g.(x) = c i-.. y i = 1 , 2,3 (B.1) 

T - D M / V 

g^(x) = u ^ 

1 - N 

where the values Ki^ M and N are the f unctioi s given by 
(B. 2) . 

K. 

1 

N 
M 




Recall that GP is the matrix shown in (B.3). 



B^i 


B '^1 


B^i 


B^i 






^-3 


^-4 


B ^2 


B ^2 


B ‘^2 


B ^2 


Bci 




BC 3 




B^3 


B-^b 


B% 


B^3 




SS 




S-4 


B^4 


B^4 


B^4 


B^4 


lifi 






aiJ 



Then given the definitions above/ the folio/ ing deriva 
tives are the result of straight forward, though tedious 
differentiation, and are offerred without detailed proof. 
LEMMA 1 



h M 



VK{T- K. )+C.T.UD 
1 11 



^ V "'i 



i = 1,2,3 <B.4) 



LEMMA 2 



LEMMA 3 






K , 
3 




V U) 



Bn 




i = 1,2,3 



LEMMA 4 



B N 
B^ 



I 



(D - V U) 



(B.5) 



(B.6) 



(B.7) 
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FCEiaUIA 1 



The first three diagonal elements of GP are the deriva- 
tives of each gi with respect to Ci (i=l,2,3) aid are given 
by (B.8) . 






. 1 



T . UDM 
1 



VK.(Ti 



V M 



2 





(B.S) 



FCRIlUlA 2 

The off diagonal elements in the first three rows and 
columns of GP are the derivatives of each gi with respect to 
Cj, and are given by (B.9). 







i = 1,2,3 
j = 1,2,3 

j f i 



(B.9) 



FCEMULA 3 

The derivative of each gi with respect to 0 is given by 
(B . 1 0 ) . 



Hi 

bi’ 



T . M ( VU - DC . ) + VK . { T . 

1 1 11 





(B. 10) 



FORMULA 4 

The first three elements of the last row of GP are the 
derivatives of g4 with respect to each Ci, and ire given by 
(B . 11 ) . 



b‘54 

be, 



9 m 

— . _ 

V (1 - N)^ 



(1-N) 



(B. 11) 
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FCBMUIA 5 

The last row, last column of 3P is the derivative of g4 
with respect to U, and is given by (B.12). 





^ ^ ''■’’'4 - (B.12) 

V (1 - N) ^ 
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APPEN^X C 

ONIFORH SAMPLES ON A TRONCATED HEMISPHERE 



The single array simulation of Chapter V required a 
random sample of positions in space uniformly distributed 
over the surface of a truncated hemisphere. Th2 hemisphere 
is tc be of radius r (3000 feet) about the acoustic center 
of a hydrophonic array. The truncation of the overall 
sphere is due to the fact that the upper portion of the 
hemisphere is above the water surface, and the lower half is 
below the sea bottom. 




Figure C.1 Hemispherical Geometry. 

Let E be the variable denoting angles of eleration above 
the horizontal. Let H be the variable denoting horizontal 
azimuth angles about the center of the hemisphere. Consider 
a small piece of hemispherical surface area bounded by the 



elevation angles El and (El+dS), and by the azimuth angles 
HI and (HI+dH) (see figure C,1). If W is the horizontal 
width of that piece, then is given by (C.1) . 

w = dH = r cos(E) dH (C.1) 

If A1 is the area of that piece, then A1 is i pproximated 
by (C.2). 



= w dE = r cos(£) dH dE (C.2) 

Now suppose that ( 0 < El < E2 < Emi x ) where 

Emax is the elevation angle of the top of the truncated 
hemishpere. Then the ratio (A1/A2) of two different areas 
at elevation angles El and E2 is given by (C. 3) . 

Aj^ r cos(E^) dH dE cos(E^) 

A^ ” r cos(E 2 ) dH dE ^ cos (E^) (C.3) 

If n1 and n2 are to be the (relative) sampls sizes from 
the two areas A1 and A2 respectively, then uniformity of the 
sample requires that (C.4) hold. 



\ / "i 



A^ / 

2 ' 2 



(C.4) 



The combination of (C.4) and (C.3) implies (1.5). 

cos (E^) 

^2 cos(Ej^) (C.5) 

Letting El =0, then the relative sample size n2.for 
area A2 is given by (C.6), where n is the relative sample 
size at the base of the hemisphere. 



= n cos (E^) 



(C.6) 



Now the differential probability of drawii g a random 
position that has elevation angle S is given by (C.7), where 
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N is the total sample size to be drawn from the surface orx 
the sector bounded by the azimuth angles HI and (Hl+dK). 






f_(E) = 

E 



COS (E) 



(C.7) 



Therefore# if K is defined to be the constant ratio n/N, 
then the corresponding cumulative distribution is given by 
(C, 8) . 

= / K cos (e) de = 



Fe(E) 



K sin (E) 



(C.8) 



For the highest elevation angle# Emax# thj cumulative 
distribution function must equal one, so that (C.9) holds. 



) 

E max 



K sin(E ) 
max 



(C.9) 



As a result# the constant K is determined by (C.10).. 

1 

(C.10) 



K 



sin(E ) 
max 



Therefore the complete cumulative distribution function 
is given by (C . 1 1) . 



Fe(E) 



sin (E) 



sin (E ) 
max 



(C. 11) 



Now the inverse probability transform can be used. If U 
is a uniform random variable on the interval (0#1)# then let 
E be given by (C.12). 



E = arcsin 



in(—^ ) 

Vsin(E )J 



max 



(C. 12) 



Now choose an azimuth angle H randomely ai d uniformly 
over the interval (0#27T). Then (X#Y#Z) given by 

X = r CO s (H) cos (E) 

X = r si n (H) cos (E) 

Z = r sin(S) 
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is the corresponding position in spherical :oordinates. 
That position is a random position drawn from a population 
of positions uniformly distributed across the sirface of a 
truncated hemisphere with maximum elevation angle Smax. 
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APPENDIX D 

SINSLE ARRAY SIMOIATICN COMPUTER PROGRaa 



THIS FORTRAN PROGRAM SELECTS A RANDOM SAMPLE OF 
SIZE M OF 3-DIMENSIONAL POSITIONS ALL AT A RANGE 
OF F. FEET FROM THE ACOUSTIC CENTER OF A HYDROPHONIC 
ARRAY. THE ARMS OF THE ARRAY ARE ASSUMED TO BE 
ALIGNED WITH THE COORDINATE SYSTEM OF THE RANGE. 

AND IS IN A POSITION GIVEN BY THE VECTOR A. 



THE SPEED 
GIVEN BY 



OF SOUND 



PROFILE IS 
V = VV(1) 



ASSUMED LINE\ R, 

+ VV (2) * DEPTH 



FOE EACH OF THE POSIIONS SELECTED, EXACT HYDROPHONE 
TIMES ARE CALCULATED USING THE SUBROUTINE TCOMP, AND 
RANDOM ERRORS ARE ADDED TO THE EXACT TIMES. THE 
ERROR DISTRIBUTION IS NORMAL, WITH MEAN ZERO AND 
STANDARD DEVIATION SDEV. 

THE TIME VALUES ARE THEN USED TO ESTIMATE AN APPARENT 
POSITION BY TWO DIFFERENT METHODS, USING APPROPRIATE 
SUBROUTINES. THE RESULTING APPARENT POSITIONS ARE 
THEN PROCESSED BY THE SUBROUTINE TRACE TO OBTAIN 
THE CORRESPONDING ACTUAL POSITION ESTIMATES. THE 
RESULTING POSITION ESTIMATES ARE THEN COMPARED TO THE 
ORIGINAL TRUE POSITION TO SEE WHICH METHOD PRODUCED 
THE MOST ACCURATE ESTIMATE. 



THIS IS DONE 
SIX METHODS. 



FOR ALL POSSIBLE PAIRS OF THE 



INTEGER M,I. J, K.N, MET H, M ET H 1 , ME TH2, N YES . N YEST , NN OISE 
DIMENSION REL( 1000) .e AZ ( 1000) -STUFF(4000) 

DOUBLE PRECISICN P(10 00.3) ,A(3) .VV {2) ,T (1000^ 

VARi 1000|^,TN^1000, 41 ,NOISE, ^RACT 



DOUBLE 

DOUBLE 

DOUBLE 

DOUBLE 

DOUBLE 



PRECISION 

PRECISION 

PRECISICN 

PRECISION 

PRECISICN 



TEST (10 0C 



I006, 3) 



DIF (1000) 



DIFT hOOo! ,AC|3),TC (l 600 ) 

_ _ A ,Pi,AZ (1000) ,EL (lUUU) 

SEED, PEST (1000,3) 



R,Z, SDEV, 



- - 

V,FK,PHI,PHIMAX, 



SET INITIAL VALUES 
AND ERROR STANDARD 

M = 1000 
R = 3000. DO 
SDEV = 3. 464D-5 
FOEMATfZX,' NOISE 
WRITE (6, if SDEV 
WRITE (6. 999) 

SEED = 931 947. ODO 
PI = 2.DO*DARCOS(0 



OF SAMPLE 
DEVIATION 



SIZE, RANGE, 



STD. 



DO) 



)EV. 



F15. 10) 



SET 

FOR 



ARRAY POSITION AND 
SPEED OF SCUND. 



LINEAR PROFILE 



ii — 

P) = 



= O.DO 
= O.DO 
= 1300. DO 



AC 

AC 






+ 



3) - 



15.0 

15.0 

15.0 
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10 

20 

11 

100 

122 

11 1 

30 1 
200 



VV (1) = 4840. 7D0 
VV [21 = 33 14. D-5 
V = VV (1 J +VV (2) *A (3) 

SET UPPER LIMIT ON ELEVATION ANGLE, IF NECil SSAEY 



GO TO 2 0 

PHIMAX = DARS IN {A (3) /R) 

GENERATE RANDOM VALUES FOR ANGLES OF 
ELEVATION AND AZIMUTH. 

FK = 1.D0/DSIN (PHIMAX) 

CALL GGUBS (SEED,M,RAZ) 

CALL GGU3S (SEED,M,REL;i 

CONVERT ANGLES TO 3-D COORDINATES 

DO 11 I = 1,M 



COS (AZ (I) *PI*2.D0 = 

SIN (AZ (I) *PI*2.D0) 

P(I,3) = A (3)-R*DS IN (PHI) 

CONTINUE 

COMPUTE EXACT HYDROPHONE TIMES UNDER THE LC NEAR 
VELOCITY PROFILE ASSUMPTION. 

CALL TCOMP (P, AC,M,VV, TC) 

CALL TCOMP (P,A,M,VV, T) 

GENERATE AND ADD ERRORS TO TIMES. 

NNOISE = 4*M 

CALL GGNML (SEED, NNOIS E, STUFF) 

K = 1 

DO 1 1 1 I = 1 ,M 
DO 122 J = 1,4 

NOISE = SIUFF(K) 

NOISE = NOISE*SDEV 
TN(I,J) = T(I,J)+NOISE 
K = K + 1 

CONTINUE 
CONTINUE 

FORMATfZX.' { SAMPLE SIZE = »,I),' )» 

WRITE (6,301) M 
WRITE (6,999) 

TUO OUTER LOOPS RUN THROUGH ALL PAIRS OF TIE SIX 
POSSIBLE POSITION ESTIMATING METHODS. 

DO 1111 METH1 =1,5 
KMETH = MET HI * 1 
DO 1222 METH2 = KMETH, 6 

IF (METHI .EQ. METH2) GO TO 1222 



FOEMAT(2X. ’METHODS USED ARE : ’) 
WRITE (6, 2 OO) 

METH = METHI 
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130 
201 

131 
202 

132 

203 

133 

204 

134 

205 

135 

206 



136 

137 



CAIL SUBEOUTINIS TO IMPLEMENT THE METHODS 

IF (METH . NE. 1) GO TO 131 

CALL POSNAV (TN , V, M, P EST , TES I) 

FOEMAT/2X.’ NAVY^ 0 NCORRE; TED’ ) 

WRITE (6-20 1 ) 

GO TO 150 

IF (METH .NE. 2) GO TO 132 

CALL POSNVC (TN,V, M, PEST, TEST) 

FORMAT (2X.* NAVY, CORRECTS D') 

WRITE (6,202) 

GO TO 150 

IF (METH .NE. 3) GO TO 133 

CALL POSLS -(TN,V,H , PEST, TEST) 

FORMAT (18X, ’LEAST SQUARES, UNCORRECTED') 
WRITE (6,203) 

GO TO 1 50 

IF (METH . NE. 4) GO TO 134 

CALL POSLSC (TN,V, M, PEST-TEST) 

FORMAT (18X, 'LEAST SQUARES, CORRECTED') 

WRITE (6,204) 

GO TO 150 

IF (METH .NE. 5) GO TO 135 

CALL POSMLP (TN,V, M, PEST, TEST, VAR) 

FOEMAT(18X. 'MAX. LIKELIHOOD, PLANAR’) 

WRITE (6,205) 

GO TO 1 50 

IF (METH . NE. 6) GO TO 136 

CALL POSMLS (TN , V. M-P EST -TEST, VAR ) 

FORMAT ( 18x- 'MAX. LIKELIHOOD, SPHERICAL' ^ 
WRITE (6,206) 

GO TO 150 



WRITE(6-137) METH1,METH2 

FORMAT(2X, 'CANT FIND METHODS *,I4,' AND/OR 
GO TO TOOO 



/I4) 



THE APPARENT POSITION ESTIMATES ARE RELATIVE TO 
THE ARRAY CENTER, AND MUST BE TRANSLATED TO 
THE TRACKING RANGE COORDINATE SYSTEM. 



150 


DO 144 I 




1,M 








PEST 1 


[l/ 


11 


= 


PEST 


(I 




PEST 


I, 


2 


1 = 


PEST 


(I 


144 


PEST I 
CONTINUE 




3 


1 ~ 


A (3) 





(1,1) + A (1 
I,2[ + A 2 
- PEST (1,3 



156 

155 



CONVERT APPARENT POSITIONS TO ACTUAL ESTIMl TES BY 
RAY TRACING PROCEDURE. 

CALL TRACE (PEST,PT, TEST,A, M, VV) 

CALCULATE DIFFERENCE BETWEEN THE 1ST POSITION 
ESTIMATE AND THE TRUE POSITION. 



IF ( METH .NE. METH1 
DO 155 I = 1,M 
DIF (I) = O.DO 



) GO TO 160 



DIFTfl) = DABS (TEST (I) -TC (1,4) ) 

DO 156 J = 1,3 

DIF(I) = DIF(I) + (PT (I,J)-P (I, J) ) **> 
CONTINUE 
CONTINUE 
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c 



160 



167 



165 

166 
C 

300 

303 

302 



1222 
11 1 1 
C 

998 

999 
1000 



MEIH = METH2 
GO TO 1 30 

CALCULATE DIFFERENCE BETWEEN THE 2ND POSITION 
ESTIMATE AND THE TRUE POSITION, AND COMPARS I^ 
TO THE 1ST POSITION DIFFERENCE. 



NYES = 0 
NYEST = 0 
DO 166 I = 
DO 167 J 



DI F f I) 
JUE 



1/M 

= f)IF(I) -(PT(I,J)-P(I,J))*=^<2 

CO NTI NU E 

DIFT(I) = DIFT(I) - DABS (TEST (I)-TC (1,4) ) 
IF'( DiFT(I) .GT. O-DO ) GO TO 165 
NYEST = NYEST + 1 
IF ( DIF(I) .GT. O.DO ) GO TO 166 
NYES = NYES + 1 
CONTINUE 



FEACT = (DFLOAI (NYES) /DFLOAT 
FORMAT (2X,» FRACTION OF TIME METHOD #1 IS 

FORMAT 1 



2X, 



FORMAT i2X. 

WRITE (6,300) 

WRITE 6,303] FEACT 
FRACT = (DFL0AT(NY3 
WRITE (6,302) FRACT 
WRITE (6, 999 
WRITE 6, 998) 

WRITE (6, 999) 



IN 

IN 



POSITION 

TIME 



: LOSER* ) 
• ,F7. 3) 

*, F7.3 



:ST)/DFLOAT (M) ) 



CONTINUE 

CONTINUE 

FORMAT (1 X, • = 
FORMAT (2X, ' 
STOP 
END 



') 



•') 
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APPENDIX E 

DOUBLE AEEAI SIMULATION COMPUTEE PBOGEIM 



1 



C 



THIS FORTEAN PROGRAM SELECTS A RANDOM SAMPLE OF 
SIZE M OF 3-DIMENSIONAL POSITIONS IN A BOX RUNNING 
CROSSWAYS BETWEEN TWO HYDROPHONIC ARRAYS. THE BOX 
HAS DIMENSIONS GIVEN BY THE VECTOR BOX. THE ARRAYS 
ARE SEPERATED BY 7500 FEET. AND HAVE COORDINATES 
GIVEN BY THE VECTORS A1 AND A2. THE ARMS OF THE 
ARRAYS ARE ALIGNED WITH THE RANGE COORDINATE SYSTEM. 

THE SPEED OF SOUND PROFILE IS ASSUMED TO BE LINEAR, 
AND IS GIVEN BY V = VV(1) + VV(2)*DEPTH . 



FOE EACH OF THE POSITIONS SELECTED, EXACT TIME VALUES 
FOE SOUND WAVS ARRIVAL AT THE HYDROPHONES OF EACH 
ARRAY ARE COMPUTED BY THE SUBROUTINE TCOMP. THEN 
RANDOM ERRORS ARE ADDED TO ALL TIMES. THE ERROR 
DISTRIBUTION IS NORMAL, WITH MEAN ZERO AND STANDARD 
DEVIATION SDEV. 



THE TIME V 
POSITIONS 
APPEOPRIAT 
DIFFERENT 
POSITIONS 
ESTIMATES 
OF THE DIF 
THE TWO DI 
METHOD PRO 



ALUES ARE THEN USED TO ESTIMATE IPPAEENT 
BY TWO DIFFERENT METHODS, USING THE 
E SUBROUTINES. EACH METHOD PRODUCES TWO 
ESTIMATES, ONE FOR EACH ARRAY. THE APPARENT 
ARE THEN TRANSLATED TO ACTUAL POSITION 
BY THE SUBROUTINE TRACE. THEN THE LENGTH 
FERENCE VECTOR FOR EACH METHOD IS COMPUTED. 
FFERENCE LENGTHS ARE COMPARED TO SEE WHICH 
DUCES POSITION PAIRS IN CLOSEST AGREEMENT. 



THIS COMPARISON IS 
THE SIX METHODS. 



DONE FOR ALL POSSIBLE Pi IRS OF 



INTEGER M ,1, J, K,METH, METH 1 , ME TH 2 , N YES , N NOI3 E,NARRAY 
DIMENSION RP (1 000) , ST OFF (4 000) , NAM 1 (6) , NAM 2 (6) 



DOUBLE PRECISION A 1(3) 
DOUBLE PRECISION VAR ( ' 
DOUBLE PRECISICN TEST 



) ,A2(3) ,T1 n000,4) ,T2{ 1 
100 0).,NOISE-FRACT,XYZ,D 
(1000) ,P1 (1000,3f,P27f 0 



DOUBLE PRECISION SDEV , SEED ^ ^OX l3)'''i^ PEs!t'( lOOO , 3) 
DOUBLE PRECISICN T (10 00, 4) , V, VV (2‘ 



(2) ,DIF(1000) 



000,4) 
DEPTH 
000,4) 
/A(3) 



DATA NAMI (1 
DATA NAM 1 (4 
DATA NAM2 (1 
DATA NAM 2 (4 



,NAM1 
,NAM1 
, NAM2 , 
,NAM2 



, NAMI 1 


[31 


1 , NAMI I 


6 


, NAM2 


3 


, NAM2i 


6 



/’NAVY’ , ’ 



NAVY’ 



L.S. 






/’L.S. ’ , 'M 

/'CORE', ’PLNR', 'SPHR'/ 



I ’ L- 

’CORE’'’ *’/ 

I ■ 



INITIALIZE VALUES FOE SAMPLE SIZE, 
AND ERROR STANDARD DEVIATION. 



RANGE, DEPTH 



M = 1000 
SDEV = 3.464D-5 
DEPTH = 1300. DO 
FORMAT f2X,’ NOISE 
WRITE (6, if SDEV 
WRITE (6,999) 

SEED = 9347.6D0 



STD. DEV. = ’,F15.10) 
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30 



22 

11 



45 

44 

C 



56 

55 



C 



SET 0? ARRAY POSITIONS, AND SOUND VELOCITY 

A1 (1) = -3750 . DO 
A1 2) = 0. DO 
A1 (3) = DEPTH 
A2 (1) = 3750. DO 
A2 (2) = 0. DO 
A2 3) = DEPTH 
VV (1) = 4840 . 7D0 

VV (2) = 33 14.D-5 
V = VVn ) +VV (2) * 



( 1 ) 

FORMAT (2X, ’ 
WRITE (6,30) M 
WRITE (6, 999) 



(2) *DEPTH 



( SAMPLE SIZE = ’,1) 



DRAW UNDERWATER BOX FOR RANDOM POSITIONS 



BOX 



BOX (2 



BOX 



= 1000. DO 
= 5000. DO 
= DEPTH 



GENERATE RANDOM POSITIONS IN BOX 

DO 1 1 J = 1,2 

CALL GGUBS (SEED,M,RP) 

DO 22 I = 1,M 
XYZ = RP (I) 

P1(I,J) = (XYZ- 0.5D0) *BOX (J) 

CONTINUE 
CONTINUE 



CALL GGUBS (SEED, M,RP) 

DO 33 I = T,M 
XYZ = RP(I) 

PI (I, 3) = XYZ*BOX(3) 

CONTINUE 

COMPUTE EXACT HYDROPHONE ARRIVAL TIMES, AND 
RANDOM ERRORS TO TIME VALUES. 

NNOISE = M*4 

CALL TCOMP (PI , A1,M,VV ,T1) 

CALL_^GGNML (SEED, NNOI SE, STUFF) 

DO 44 I = 1,M 

DO 45 J = 1,4 

NOISE = STUFF(K) 

NOISE = NOISE*SDEV 
T1 (I-J) =T1(I,J) + NOISE 

K = K+1 
CONTINUE 
CONTINUE 

CALL TCOMP (PI , A2, M,VV ,T2) 

CALL GGNML (SEED, NNOI SE, STUFF) 

K = 1 

DO 55 I = 1,M 
DO 56 J = 1,4 

NOISE = STUFF(K) 

NOISE = NCISE*SDEV 
T2 (I,J) = T2 (I, J) + NOISE 

K = K+1 
CONTINUE 
CONTINUE 



PROFILE 






ADD 



100 



c 

c 

c 

10 

20 

C 

c 

c 

c 

125 

67 

66 

C 

C 

130 

131 

132 

133 

134 

135 

136 

137 

C 

150 

C 

C 

C 

C 

c 

144 

C 

C 

c 

c 



SET UP LOOPS TO EUN THHOUGH ALL PAIRS OF MSTHODS. 

DO 1111 METH1 =1,5 
KHETH = HETH1 + 1 
DO 1222 METH2 = KMETH.S 

IF (METH1 .EQ. METH2) GO TO 1222 

FORMAT {2X, ’METHODS COMPARED ARE 1, ’,A4,2X,A4) 

F0RMAT]2X,’ and 2. ’,A4,2X,A4) 

WRITE (6, 10) NAM1 (METH 1) ,NAM2(METH1) 

WRITE (6, 20[ NAM 1 (METH 2) ,NAM2(METH2 
WRITE (6,9 9 9) 

METH = METH1 
NAERAY = 1 

SET DP ARRAY BEING USED 



^ ( 


1) 


= A1 1 


■I) 


A 


2) 


= A1 1 


2 


A 


3) 


= A 1 1 


3) 


DO 


66 


-L — 


1,M 



CONTINUE 

CALL SUBROUTINES TO PERFORM ESTIMATION MET30DS. 

IF (METH .NE. 1) GO TO 131 

CALL POSNAV (T , V, M , PEST , TESI ) 

GO TO 1 50 

IF (METH . NE. 2) GO 1 0 132 

CAIL POSNVC (T,V,M , PEST, TEST) 

GO TO 150 

IF (METH .NE. 3) GO TO 133 

CALL POSLS (I, V,M, PEST, TEST) 

GO TO 1 50 

IF (METH .NE. 4) GO TO 134 

CALL POSLSC (T, V,M , PEST, TEST) 

GO TO 1 50 

IF (METH ,HE. 5) GO TO 135 

CALL POSMLP (T,V,M , PEST, TESI, VAR) 

GO TO 150 

IF (METH .NE. 6) GO TO 136 

CALL POSMLS (T, V, M , PEST, TEST , VAR) 

GO TO 150 

WRITE (6. 137) METH1.MEIH2 

F0EMAT(2X, ’CANT FIND METHODS ’,14,’ AND/OR ’,14) 

GO TO 1000 

IF ( NARRAY . EQ. 2 ) GO TO 152 

APPARENT POSITION ESTIMATES ARE IN LOCAL AJRAY 

COORDINATES, AND MUST BE TRANSLATED TO TRACKING 

RANGE SYSTEM COORDINATES. 

DO T = 1 M 



CON 

CORRECT APPARENT POSITIONS BY RAY TRACING. 
CALL TRACE (P EST, P 1, T EST ,A , M, VV) 
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GO ON TO SECONE ARRAY 



156 

155 



(1) 


= A2 ( 1) 


2) 


= A2 2) 




^2 3 



NARRAY = 2 
A 
A 
A 

r, H 

DO 78 J = 1 , 4 

T(I, J) = 12 (I, J) 

78 CONTINUE 

77 CONTINUE 
GO TO 13 0 

TRANSLATE 2ND ARRAY APPARENT POSITIONS TO I ANGE 
COORDINATE SYSTEH, AND THEN RAY TRACE. 



152 


DO 145 I 


: = i,M 






PESTj 


[1,1) = 


PEST (I 




PEST 


;i,2) = 


PEST (I 




PESTi 


1,3 = 


A2(3) 


145 


CONTINUE 



,1) + A2 (1) 

,2) + A2(2 
- PEST (i;s 



CAIL TRACE (PEST, P2, T EST , A , M, VV ) 

COMPUTE DIFFERENCE 7ECTORS FOR 1ST METHOD 

IF { MEIH .NE. METH1 ) GO TO 160 
DO 155 I = 1,M 
DIF (I) = O.DO 
DO 156 J = 1,3 

DIF(I) = 5iF(I) + (P1 (I,J)-P2(I,J))**2 
CONTINUE 
CONTINUE 

GO ON TO SECOND METHOD 

METK = METH2 
NARRAY = 1 
GO TO 125 



COMPUTE DIFFERENCE VECTORS 
COMPARE TO DIFFERENCES FOR 



FOR 2ND METHOD i ND 
1ST METHOD. 



160 


NYES = 


: 0 




DO 166 I = 




DO 


167 J 
DIF (I) 


167 


CONTINUE 




IF 


( DIF( 
NYES = 


,166 


CONTINUE 




WRITE (6,999) 




FRACT 


= (DFL 


300 


FOEMA"] 


: (2 X , ’ 


310 


FORMA1 


:j2x'’ 
[6, 300) 




WRITE 




WRITE 


6,310) 




WRITE 1 


6, 99 9) 




WRITE 1 


6,998 




WRITE ( 


6, 99 9) 



= 1,3 



.GT. 

I YES + 



0 

1 



DO ) GO TO 166 



C 

1222 

11 1 1 

C 

99 8 
999 
1000 



OAT(NYES) /DFLOAT(M) ) 
FRACTION OF TIME METHOD 



IS BETTER IS 



#1 M 
• ,F7.3- r 



FRACT 



CONTINUE 

CONTINUE 

FORMAT (1 X, 
FORMAT (2X, ' 
STOP 



’) 
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AP PEN DIX F 

COMPOTEE SOBROailNES FOR TIME CALCOLATION AID RAITRACING 



SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS3 SSSSSSSS 
SUBROUTINE TCOMP (P,A,M,VV,T) 



22 

11 



44 

33 



THIS FORTRAN SUEEOUIINE COMPUTES THE EXACT TIME 
REQUIRED FOR A SOUND WAVE TO TRAVEL FROM ITS 
SOURCE (VECTOR F) TO THE FOUR HXDROPHONES 0N AN 
ARRAY WHOSE ACOUSTIC CENTER IS SPECIFIED BY 
VECTOR A. 



THE METHODOLOGY USED IS SET FOETH IN APPENDIX A 
OF THE THESIS. THE BASIC ASSUMPTION IS ONE OF 
A LINEAR VELOCITY PROFILE, WHOSE COEFFICIENTS ARE 
GIVEN BY THE VECTOR VV. 

INTEGER M,I,J 

DOUBLE PEECISICN P ( 1 000,3) , A (3) , VV ( 2) . T ( ID 0 0, 4) 
DOUBLE PRECISION AA ( 4, 3) , P 1 , P2 , R, C 1 , A2, C2 



SET UP COORDINATES OF FOUR HYDROPHONES ON THE ARRAY 
DO 

DO 22 J = 1.3 

DO + A (J) 



111= 1,4 
DO 22 J = 1,3 
AA (I,J) = -15, 

CONTINUE 

AA^I^3) = 15. DO * A (3) 



+ A(3) 



CONTI 

AA (1,1) = 15. DO 
AA 12,2) = 15. DO 
A A (3,3) = -15. DO + A (3) 

CALCULATE QUANTITIES Cl, C2, R, AND THE TI! ES T(*,4) 

C2 = -1.D0*(VV (1)/VV(2) ) 

lOCP THROUGH THE M SOURCE LOCATIONS 
DO 33 I = 1,M 
P2 = P (1,3) 

LOOP THROUGH THE FOUR HYDROPHONES FOR Ek CH SOURCE 
DO 44 J = 1,4 

PI = DSQRT ( (P(I ,11-AA (J, 1) ) **2 

+ (f (l,2)-Ai(J,2) ) **2) 

Cl = P1**2+P2**2-AA (J, 3) 

+ 2. DO *C2* (AA(J, 3) -P2) 

Cl = Cl/ (2.D0 + P 1) 

RT ' 



R = DSQR' 
T(I, J) 



CO NT 

CONTINUE 



T(I, J) 
INP- 



UE 



(C1*^2+ (AA ( J,3) - C2) **2) 
((AA(J,3)-C2)/(P2-C2) ) 

* ( (R-C1+P1 / (R-CI ) ) 
(DLOG(T(I,J)) )/VV(2) 



RETURN 

END 
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SS5SS5SSSSSSSSSSSSSSSSSSSSSSS3SSSSSSSSS3SSSS3SSSSSSSSSSS 
SUBROUTINE TRACE (P , PT , T , A, M, V V) 

THIS FORTRAN SUBROUTINE TAKES AN APPARENT POSITION 
ESTIMATE P AND CONVERTS IT TO AN ACTUAL POSITION 
ESTIMATE PT BY RAY TRACING. THE RAY TRACING ASSUMES 
A LINEAR VELOCITY PROFILE, WITH COEFFICIENTS IN THE 
VECTOR VV. OTHER INPUTS ARE THE ARRAY LOCATION 
VECTOR A AND THE RAY TRANSIT TIME TO THE ACOUSTIC 
CENTER, T. THE RAY TRACING LAYERS ARE DEL FEET THICK 
WITH THE FIRST (BOTTOM) LAYER BEING THE FIP ST DEL 
FEET IMMEDIATELY ABOVE THE ARRAY. 



10 



INTEGER H, I 

DOUBLE PRECIS ICN T (10 00) ,P( 10 00,3) ,PT (1000, 3) ,A (3) 
DOUBLE PRECISION DEL , H ,S A, C A, DEPTH , V , V 1 , V V ( 2) 
DOUBLE PRECISION DS,DT,DZ,DH 

DEL = 25. DO 

LOOP THROUGH THE M APPARENT POSITIONS 
DO 11 I = 1,M 

INITIALIZE INCREMENTAL VALUES FOR EACH POSITION 
H = DSQRTf (P (I, 1)^- A (1)^) **2+ (P (1,2) -A (2) I **2) 

CA = H/DS0RT(H**2+ (P (I,3)-A(3) ) **2) 

DT = O.DO 
DZ = 0. DO 
DH = 0. DO 

DEPTH = A (3) - DEL/2. DO 

V = VV (1) +VV (2) *DEPTH 

INNER LOOP : PROCESS DATA UPWARDS, LAYER BY LAYER, 
UNTIL RAY TRANSIT TIME IS EXHAUSTED 

SA = DSQRT f l.D0-CA**2) 

DZ = DZ DEI 
DS = DEL/SA 
DT = DT DS/V 
DH = DH + DS*CA 

IF ( DT .GT. T (I) ) GO TO 20 

DEPTH = DEPTH - DEL 
VI = VV(1) + VV(2) *DEPTH 

CA = CA*(V1/V) 

V = VI 
GO TO 10 

ADJUST FOR OVERSHOOTING IN LAST LAYER 



20 DS = V* (DT-T (I) ) 

DH = DK - CA*DS 

ADJUST APPARENT POSITION TO GET ACTUAL POSC TION 



11 



PT (I, 1 
PT (I, 2 
PT (I, 3 

CONTINUE 

RETURN 

END 



= ] 


?(I,1)-A(1) 


) *(DH/H) 


+ A ( 


= 


?]l,2 -AJ2 


) *(DH/H 


+ A ( 


= A 


(3) - (DZ-S 


A*DS) 
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APPENDIX G 

COMPUTES SUBROUTINES FOR POSITION ESTIMATION METHODS 



22 

T1 

C 



;SSSSSSSSSSSSSSS£SSSSS SSSSSSSSSSSSS3SSSS SSSSSSSSS3SSSSS 
3UBR0UTINE P03NAV {T,V,M,P,NT) 

INTEGER M, I 

THI3 FORTRAN 3U3ROUTINE IMPLEMENT3 THE ORI5 INAL 
UNADJU3TED NAVY APPARENT P03ITI0N E3TIMATI0N METHOD. 

INPUT ARE THE HYDROPHONE TIMES T FOR M SOU^ D SOURCE 
POSITIONS, AND VELOCITY V AT THE ARRAY. 

OUTPUT ARE M APPARENT POSITION ESTIMATES P ALONG WITH 
THE CORRESPONDING M RAY TRANSIT TIMES NT. ALL THE 
POSITIONS ARE REFERENCED TO THE ACOUSTIC C2NTER, WITH 
THE Z COMPONENT BEING MEASURED UPWARDS FROM THE ARRAY. 

DOUBLE PRECISION T ( 10 00 , 4) , P { 1 0 00 , 3) , NT (1 0) 0) 

DOUBLE PRECISICN V , TC , D, BIS C, NU MER, D ENOM 

D = 30. DO 
DO 11 I = 1.M 
NUMER = O.DO 
DENOM = O.DO 



CONTINUE 

NT (I) = TC*DSQRT(NUMER/DEHOM) 

CONTINUE 
FETURN 
END 
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SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSj sssssssssss 
SUEEOUTINE POSNVC {T,V,M,P,NI) 

THIS FORTRAN SDBROUTINE IMPLEMENTS THE ADJUSTED NAVY 
METHOD NAVY_A FOR ESTIMATING APPARENT POSITIONS. 

INPUT ARE THE HYDROPHONES TIMES I FOR M SOJ ND SOURCE 
LOCATIONS, AND THE VELOCITY V AT THE ARRAY. 

OUTPUT ARE THE H APPARENT POSITION ESTIMATES P, AND 
THE CORRESPONDING M RAY TRANSIT TIMES NT. ALL 
POSITIONS ARE REFERENCED TO THE ACOUSTIC CENTER OF 
THE ARRAY, WITH THE Z COMPONENT MEASURE UPWARDS 
FROM THE ARRAY. 



22 



44 

11 



INTEGER M,I 
DOUBLE PRECISICN 
DOUBLE PRECISICN 



T(1 000,4) ,P (1000,3) , NT (13 00) 
V,D,NUMER,DENOM, TC,DCC 



D = 30. DO 
DO 1 1 I = 1, M 
DCC = O.DO 
TC = T(I,4) 

DO 22 J = 1,3 

P(I,J) = 15. DO + {V*V* (TC*TC-T {I, J) *1' 2) / (D*2. CO) ) 
DCC = DCC + P(I,J)**2 
CONTINUE 
NUMER = O.DO 
DENOM = O.DO 
DO 44 J = 1,3 

PJIcJ) = (P (I, J) *V*TC)/DSQRT (DCC) 

DENOM = DENOM * P(I,J)**2 
P(I,J) = P(I,J) - 15. DO 
NUMER = NUMER P{I,J)**2 
CONTINUE 

NT (I| = IC*DSQRT(NUMER/DENOM) 



CONTIN 

RETURN 

END 
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ssssssssssssssssssssssssssssssssssssssssssssssssssssssssss 

SUBROUTINE POSTS { T, V , M , P, LSI) 

THIS FORTRAN SUBROUTINE IMPLEMENTS THE LEAST SQUARES 
PLANAR METHOD L.S. 

INPUT ARE THE HYDROPHONE TIMES T FOR M SOUND SOURCE 
POSITIONS, AND THE VELOCITY V AT THE ARRAY. 

INPUTS AND OUTPUTS ARE THE SAME AS FOR THE 
SUBROUTINE POSHAV. 

INTEGER M-I-J 

DOUBLE PRECISION T ( 1 0 00 , 4) , p ( 1 0 00 , 3) , LST ( 1 3 00 ) 

DOUBLE PRECISICN V,DISC,TC 

DO 1 1 I = 1, M 
DISC = O.DO 
DO 22 J = 1,3 

P(I,J) = T(I/4) -T(I,J) 

DISC = DISC + P(I,J)**2 



(V*LST ( I) *P (I, J) ) /DSQRT (DISC) 



22 


CONTINUE 




LST(I) = 
DO 33 J = 


33 


P(I.J) 

CONTINUE 


11 


CONTINUE 




RETURN 

END 
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SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS33SSSSSSSS5 sssssssssss 
SUBROUTINE POSISC (T , H , IS T) 

THIS FORTRAN SUBROUTINE IMPLEMENTS THE BIAS ADJUSTED 
LEAST SQUARES PLANAR METHOD L.S.C. FOR ESTIMATING 
APPARENT POSITIONS. 



55 

44 



22 



33 



89 



88 



77 



66 



99 



INPUTS AND OUTPUTS ARE THE 
SUBROUTINE POSTS. 



SAME AS FOR THE 



INTEGER I.M.J.K.L 

DOUBLE PRECISICN 1 (10 00.4) ,P( 1000^3) ,LST(1) 00) 

DOUBLE PRECISION V , A ( 4, 3 ) , D (4 ) , DISC, TD, R IS QR , R2SQR, Z (3) 

SET UP HYDROPHONE POSITIONS 
DO 44 J =1,3 

DO 55 I = 1,4 

A(I-J) = -15. DO 
CONTINUE 
AiJ,J) = 15. DO 
CCNTINDE 

CALCULATE ORIGINAL L. S. SOLUTION 
DO 11 I = 1,M 
DISC = O.DO 
DO 22 J = 1,3 

F(I,J) = T(I,4) -T(I.J) 

DISC = DISC ♦ P(I,J)**2 
CONTINUE 

LSTil) = (T(I,1)+T(I,2)+T(I,3)-T(I,4))/2. 

DO 33 J — 1 # 

= (V*LST(I) *P (I, J)) /DSQRT (DISC) 

CONTINUE 

R1SQE = P(I, 1)**2+P (1,2) **2 + P (1,3) **2 
D (4) = O.DO 
DO 89 K = 1,3 

^^„gi4) = D(4) + (P(I^K) - A(4,K))=!<*2 

CONTINUE 

D (4) = DSQRT (D(4) ) 

DISC = O.DO 
TD = O.DO 
DO 77 J = 1 , 3 
D(J^ = O.DO 



DO 



DO 88 K = 1,3 
D(J) = D(J) 
CONTINUE 



+ (P (I,K) - A (J,K)) **2 



D(J) = DSQRT (D(J) ) 

DISC = DISC + (D (4) -D (J) ) **2 
TD = TD + D(J) 

CONTINUE 

CALCULATE BIAS VECTOR E 
DISC = DSQRT (DISC) 

DO 66 J = 1,3 

= ((TD-D(4) )*(D(4)-D(J))/(DISC*2. ,D0))-P(I,J) 

CCNIIND 
R2SQR = 



DO 



0. DO 
9 J =1,3 



R2SQR = 
CONTINUE 



R2‘ 



P(lij| 



**2 



11 



ADJUST ORIGINAL RAY TRANSIT TIME 
1ST (I) = LST (I) * (DS QRT (R2SQR/R1SQR) ) 

CONTINUE 
RETURN 
END 
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SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS3SSSSSSSSSSSSSSSSSSSSS 



12 

22 



33 



66 



60 

61 



SUBROUTINE PO SMLP (T , V , M, P, M T, VAR) 

THIS FORTRAN SUBROUTINE IMPLEMENTS THE MAXC MUM 
LIKELIHOOD PLANAR METHOD M.L.P. FOR ESTIMATING 
APPARENT POSITIONS. 

INPUTS ARE THE HYDROPHONE TIMES T FOR M SOURCE 
POSITIONS, AND THE VELOCITY OF SOUND V AT THE ARRAY. 

OUTPUTS ARE THE APPARENT POSITION ESTIMATES P AND 
RAY TRANSIT TIMES MT FOE EACH OF THE M POSITIONS. 
ALSO OUTPUT IS A VECTOR VAR OF M ESTIMATES OF THE 
TIMING ERROR VARIANCE. POSITIONS ARE ALL REFERENCED 
TO THE ACOUSTIC CENTER OF THE ARRAY, WITH THE Z 
COMPONENT MEASURED UPWARDS FROM THE ARRAY. 

INTEGER M, I, J 

DOUBLE PRECISICN T (1 0 00 , 4) , P ( 10 00 , 3) , MT ( 1 03 0) 

DOUBLE PRECISION V, TC , D , DIFF, NU MEfi ,D ENO , TO L , C (3) 
DOUBLE PRECISICN VAR ( 1 00 0) , CO (3) , X (3 ) 

D = 30. DO 
TOL = 1.D-4 

INITIALIZE THE DIRECTION COSINE ESTIMATE Bf 
USING THE LEAST SQUARES SOLUTION. 

LOOP THROUGH THE M POSITIONS 

DO 11 I = 1,M 
TC = T (1,4) 

DENOM = (TC-T (1,1) ) *=<'2+ (TC-T (I, 2) ) **2 

+ (TC-T I, 3) **2 

DO 12 J = 1,3 

C(J) = (TC - I (I, J)) /DSQRT (DENOM) 

CONTINUE 

ITER = 0 
DENOM = O.DO 
ITER = ITER + 1 
TC = 0. DO 

USE ITERATION FORMULAE TO DEVELOPE COSINES 
AND TIME VALUES. 

DO 33 J = 1 ,3 
CO (J) = C(J) 

DENOM = DENOM + T(I,J)*C(J) 

TC = TC + T (I,J) 

CONTINUE 

TC = ( (TC + T (1,4) ) ♦ D* (C(1) 4-C(2) 

DENOM = DENOM - TC * (C (1 ) +C (2 ) +C (5) 

DO 66 J = 1,3 

lENOM 

CONT"”*'- ^ 



66 J = 1,3 

C(J) = (T (I, Jli-TC) /DEN 
X]j) = DABS (CO ( J)-C(J) 
JITINUE 



+ C(3j)/V) /4.D0 



CHECK ITERATION TOLERANCES 

DIFF = DMAX1 (X(1) .X(2) ,X (3)) 

IF 7 ITER .LT. 100) GO TO 70 

FORMAT (2X, 'ITERATIONS EXCEED ',16) 
FORMAT (2X, 'IN POSITION NO. '/I6) 
WRITE(6,60) ITER 
WRITE(6,61) I 
GO TO 11 
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70 



IF { TOL .LI. DIFF ) GO TO 22 



44 



55 

11 



NUMER = O.DO 
DENOa = O.DO 

CONVERT COSINES TO POSITIONS, AND TRANS. ATE TO 
ACOUSTIC CENTER REFERENCED COORDINATES. 

DO 44 J = 1 , 3 

P(I- J) = V*C (J) *TC 
DENOH = DENOM * P{I,J)**2 
P(I,J) = P(I,J) -15. DO 
NUMER = NDMER + P(I,J)**2 
CONTINUE 

MT (I) = TC*DSQRT (NUHER/DENOM) 

ESTIMATE VARIANCE 

VAR (I) = 0. DO 

DO 55 J = 1,3 

VAR (I) = VAR (I) + (T (I, J) -TC+D*C{J) /M **2 

CONTINUE 

VARjl) = VAR (I)/4. DO 



CONTINUE 
fEIURN 
END 
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SUEROUTINE POSMLS (T IMS , V , M, P, MLT , V ARHL) 

THIS FORTRAN SUBROOTINE IMPLEMENTS THE MAXC MUM 
LIKELIHOOD SPHERICAL METHOD M.L.S. FOR ESTIMATIN: 
APPARENT POSITIONS. 



INPUTS AND OUTPUTS ARE THE 
SUEROUTINE POSMIP. 



SAME AS FOR THE 



INTEGER H, I, J, ISR, II. ITERl. ITER2 




DOUBLE 

DOUBLE 

DOUBLE 

DOUBLE 

DOUBLE 



PRECISION 

PRECISICN 

PRECISION 

PRECISICN 

PRECISICN 



ML(1000 
DSDC (3 
F 



3 ,5:c,ta6,tk (5) (5) ,c {3) ,ci p) 

j3 ,DENOM,NUMER,DLDC(3) .0LDT,DSDT 
P(4,41 ,GPI (4,4) , WK ARE A (2 8) ,GC(4) 



. - ] /GPI (4,4) , ^ , 

GpXI(lt) ,X2{4) ,M4) ,B(4) ,T1 
F1,F2 



SET GLOBAL VALUES 



11 1 



R = (3.D0-DSQRT (5. DO) ) /2.D0 
D = 30. DO 
TOL1 = 1 . D-5 
TOI2 = 1.D-5 

START OUTERMOST LOOP (ONE FOR EACH SOURCE POSITION) 

DO 5 II = 1,M 
DENOM = 0. DO 
TAD = TIMEiII,4) 

EPS = 1 . D-3 
TC = TAU 
DO 111 J = 1,3 
T 
C 
D 

CONTINUE 



11 J = 1,3 

(J) = TIME(II,J)^ 

jj) = ( ( V**2) * (TC**2-T (J) 

ENOM = DENOM + C (J) **2 



** 2 ) / (D*2. D) )) +D/2.D0 



122 



10 



DO 122 
CO NT] 



J = 1,3 



(J) = C (5)/DSQRT (DENOM) 

INUE 

START MIDDLE LOOP (CREATE DERIVATIVE MATRIX GP ( ,) ) 



310 
31 1 



300 



ITERl = 0 
II = TAD 

ITERl = ITERl + 1 

IF (ITERl .LT. 50 ) GO TO 300 

FORMAT(2X, ’EXCESSIVE NEWTON ITERS.) 
FORMITOY 1 IN POSITION ’ 



FORMAIf2X, 
WRITE'(b,3 10) 
WRITE (6,3 11) 
GO TO 5 
L = 0. DO 
S = 0. DO 
DLDT = O.DO 
DSDT = O.DO 



,16) 



II 



DO 11 I = 1,3 

=~TAiS*^*2 + 

TK]I) = T (I)'-DSQRT (K 
DLDC(I) = ( (V*K (I) *TK (I) ) +C(I) *T (I) *0*“ TAU)/VK(I) 

DSDC (I) = ]T (I) *D*TAU) /VK (I) 

DLDT = DLDT+ (C(I) *T (I) * (D*C(I) -V*TAU) ) /VK (I) 



(D/V) **2-f2. D0*D*TAU*C(C ) ) /V 
V*K(I)*DSQRT(K(I)) 
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11 



33 



DSDT = DSDI+T (I)* (D*C (I) -V*TAU) /VK (I) 

L= L + (C (I)*TK(I) )/DSQSI (K (I) ) 

S = S+TK (I)/DSQRT (K (I) ) 

CONTINUE 

DO 22 I = 1,3 
DO 3 3 J = 1,3 

I,J)= (TK(I) *DLDC (J) ) /(L*L*DSQRT( K (I) ) ) 



GP 
CONTI 
GP (4 , I) =DSDC 
GP(4, I) = (GP (4' 



(V*TC-D*L) 



D*DLDC (I) * (1 .DO-S) ) 



/(V*(1. DO-S 

GP(I, I) = (1*T (I) * D*TAU) -DLDC(I) *K (I) •“ 

GPi I,I) = 1.D0 - (GP (1,1) / (VK (I) *L*L) I 
GP(I,4) = (T (I) (V*TAU-D*C U) ) |^ 



* 2 ) 



*TK (I) 

*K (I) *C K (I) *DLDT 



22 



40 

44 



GP(I,4) = GP(I,4) / (VK(I)*L*L) 

CONTINUE 

GP (4,4) =1.D0-((DSDT*(V*TC-D*L))-D*DLDT*f1.D0-S) ) 

/(V* (i DO-S) **2) 

INVEST DERIVATIVE MATRIX 

CALL GAUSS3 (4,0,GP, GPI,IER, 4) 

CALCULATE INITIAL NEWTON SEARCH SOLUTION 
DO 44 I = 1,3 

= C (I)-TK(I)/(DSQRT(K(I)) *L) 

CONTINUE 

GC(4) = TAD- (V*TC-D*L) /(V<= (1 .DO-S) ) 

B(4^ = TAU 



56 

55 



66 



DO 



1 = 1,3 



B 

B 

DO 



= B -GPI (4 ,1) *GC (I) 



IH : B®- 

5 56 J = 1, 



B(I) = B (I) -GPI (I,J) *GC(J) 
CONTINUE 
CONTINUE 

B(4) = B (4) - GPI (4 ,4) ^GC (4) 

PREPARE FOR GOLDEN SECTION SEARCH 

EPS = DMAX1 ( (1.D-6) , (EPS/1. D1 ) ) 
DO 66 I = 1,3 






R*(B(I)-A(D) 



XT 

Ci 

CONI 
A(4) = TAD 

XT (4) = A (4) +R* (B (4 ) -A (4) ) 



70 

410 

400 

77 



TAU = XI (4) 

START INNERMOST LOOP ( COMMENCE GOLDEN SECTi! ON SEARCH 
TO IMPROVE THE INITIAL NEWION SOLUTION ) 

ITER2 = 0 

CALL OBJFCN (F 1 , L, K , TK ,S , XI , D , V, T, TC) 

ITER2 = ITER2 + 1 

IF (ITER 2 .LT. 50 ) GO TO 4 00 

FORMAT(2X, ’EXCESSIVE G.S.S. ITERS., P) S #’,I6) 
WRITE (6,4 10) II 
GO TO 5 
DIFF = O.DO 
DO 77 I = 1,3 



DIFF = DiFF + (A (I) -3 (I) ) **2 

|ij -'xr(i) 



X2(I) = A 
C(I) = X2 
CONTINUE 



* 



B 
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88 

90 

99 

100 
144 



155 

5 

C 

C900 

C 



X2 (4) = A(4) + B (4) - XI (4) 

lAO = X2(4) 

DIFI = DS2RT(DIFF + (A ( 4) -B ( 4) ) **2) 

CHECK TOLEHANCBS - STOP G.S.S. IF TIGHT ENOJ GH 
IF ( EPS .GT. DIFF ) GO TO 100 
CAll 03JFCN (F2,L,K, TK^ S, X2 , D, V, T , TC) 

CHOOSE IMPROVING DIRECTION 



IF 

DO 



FI .LT. F2 ) GO TO 90 
81 = 1,3 
= XI 



A(I) = XI (I) 
Xl7l) = X2(I) 
cd) = XI (i) 

CONTINUE 
A(41 = XI (4) 

XI (4) = X2 (4) 

GO TO 70 



DO 99 I 



= 

Bfll = X2 (I 
X1(I) = A (I 

C(I) = XI (I 
CONTINUE 



+ B(I)-X1(I) 



B]4) = X2(4) 

+B(4)-X1(4) 

GO TO 70 

CHECK TOLERANCES FOR NEWTON SEARCH ITERATI) NS, 

AND PREP FOR NEXT NEWTON ITERATION IF NECESSARY 

DO 144 I = 1,3 



X(I) = DA^S (C (I) -Cl (I) ) 
TINUE 



CONT 

DIFF = DHAX1 (X (1) ,X (2) ,X (3) ) 

IF ( TOL1 .LT. DIFF ) GO TO 10 

DIFF = DABS (TAU-T11 

IF { TOL2 .LT. DIFF ) GO TO 10 

DONE WITH II-TH SET OF TIMES AND POSITIONS. MAKE 
ESTIMATES, AND GO ON TO NEXT POSITION TO BE ESTIMATED 

VARML(II) = O.DO 
DENOM = O.DO 
HOMER = 0. DO 
DO 155 J = 1,3 

VARML (II) = VARHL (II) 

P(II,J) = V*TAU*C 
DENOM = DEHOM + P 



I 



TK (J) *TK (J) 



J) **2 

=_y*?AU»C (Jj_-D/2.D0 



II, J) **2 



NUMeS = NDMER + P 
CONTINUE 

VARHL (II) = (VARML ( II) + (TC-TAO) **2) /4 .DO 
SDEV(II) = DSCRT (VA RMI (II) ) 

HLT (II) = TAC*DSQRT (NUMER/DENOM) 

CONTINUE 

WRITE(7,900) (SDEV(I),I = 1,M) 
format (5 El 5. 5) 

RETURN 

END 
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3UEE0UTINE OBJICN (F, L, K K ,3 , X / D , V , T, TC) 

■IHI3 FORTRAN 3U5E0UTINE 13 CALLED BY THE 3JBROUTINE 
P03MLS. IT CALCULATE3 NEW VALUES FOR SEVERAL 
VARIABLES AS WELL AS A FUNCTIONAL VALUE WHICH IS THE 
DECISION FACTOR DETERMINING THE APPROPRIATE 
IMPROVING DIRECTION FOE THE GOLDEN SECTION SEARCH. 

INTEGER I 

DOUBLE PRECISION F, L, K (3) , TK ( 3) , 3 , X ( 4) , D , V, TC , T ( 3) 

3 = O.DO 
L = O.DO 
DO 111= 1,3 
K(I) = 



11 



22 



■,-) = X (4) **2+ (D/V) **2- 
TK (I) = T (I)-DSQRT (K (I) ) 

~ + TK (I) /DSQRT (K (I) 
+ (X (I) *TK(I) )/DSC 



** 2 - ( (2. D0*D*X (4) *X( 1} ) /V) 



))/DSQRT (K (I) ) 



3=3 
L = L 
CONTINUE 
F = O.DO 
DO 22 I = 1,3 

F = F+ (X (I) -TK (I) / (DSQRT (K (I) ) *L) ) **2 
CONTINUE 

F = F+(X (4)- (V*TC-D*L)/(V* (1.D0-S) ) ) **2 

RETURN 

END 
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